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Chapter  1 
Introduction 


The  prediction  of  noise  generated  by  hydrodynamic  sources  is  a  significant  problem  of  interest 
to  the  U.S.  Navy.  This  proposal  describes  a  program  of  research  aimed  at  developing  and 
testing  computational  techniques  for  formulating  and  solving  structural  acoustics  problems 
driven  by  hydrodynamic  sources. 

The  typical  parameters  associated  with  the  problems  of  interest  to  the  Navy  are  low 
Mach  numbers  and  high  Reynolds  numbers.  This  regime  naturally  leads  to  the  application 
of  Lighthill’s  acoustic  analogy  to  solve  the  problem  [1,  2],  In  this  approach  the  original 
problem  is  decomposed  into  two  parts,  one  that  involves  the  solution  of  the  Navier-Stokes 
equations  to  determine  the  unsteady  fluid  variables,  and  another  that  involves  the  solution 
of  an  acoustic  problem  driven  by  quadrupole  sources  whose  distribution  is  determined  by  the 
components  of  Lighthill’s  turbulence  tensor.  This  tensor  is  constructed  from  the  unsteady 
flow  field  computed  in  the  fluid  calculation.  In  our  previous  work  [3,4],  we  have  developed 
a  methodology  that  uses  large  eddy  simulation  (LES)  to  calculate  Lighthill’s  tensor  and 
a  variational  formulation  of  Lighthill’s  acoustic  analogy  to  solve  the  acoustic  problem. 
The  effectiveness  of  this  method  was  demonstrated  in  [4],  where  we  calculated  the  noise 
generated  by  turbulent  flow  over  an  airfoil,  while  fully  accounting  for  its  geometry.  Previous 
studies  [5-11]  have  had  to  make  simplifications  to  render  the  problem  tractable.  In  the 
same  paper,  it  was  found  that  the  computational  costs  of  the  overall  methodology  were 
dominated  by  the  turbulent  calculation.  This  meant  that  better  and  more  efficient  LES 
models  were  required  to  extend  the  applicability  of  the  proposed  methodology  to  higher 
Reynolds  number.  With  this  in  mind  a  parallel  effort  aimed  at  developing  new  LES  models 
based  on  the  variational  multiscale  formulation  was  undertaken  (see  [12, 13]).  The  proposed 
research  extends  these  ideas  in  new  directions. 

A  significant  component  of  the  research  is  aimed  at  improving  the  variational  multiscale 
formulation  of  LES  [12-14].  This  is  accomplished  by  developing  a  dynamic  version  of  this 
method,  wherein  the  eddy  viscosity  is  calculated  during  the  numerical  simulation  and  is  not 
fixed  a-priori.  For  this  purpose  the  dynamic  Smagorinsky  model  [15]  is  revisited,  and  a 
variational  version  of  the  Germano  identity  is  derived.  This  is  done  in  the  context  of  an 
abstract  variational  problem  in  Chapter  2.  It  is  found  that  this  leads  to  a  formulation  which 
may  be  used  to  determine  unknown  parameters  in  a  generic  numerical  method. 

Thereafter  the  variational  Germano  identity  is  applied  to  solve  problems  posed  as 
conservation  laws.  The  specialization  of  this  identity  to  account  for  spectral  discretizations 
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is  developed  and  the  effectiveness  of  this  formulation  in  solving  problems  with  discontinuities 
(shocks)  is  demonstrated.  These  developments  are  presented  in  Chapter  3. 

In  Chapter  4,  the  variational  Germano  identity  is  applied  to  the  large  eddy  simulation  of 
incompressible  turbulent  flows,  and  it  is  demonstrated  that  it  is  more  accurate  and  robust 
than  the  filtered  counterpart. 

Finally  in  Chapter  5,  the  performance  of  this  identity  in  predicting  the  far-field  noise 
generated  by  homogeneous  turbulence  is  assessed.  It  is  found  that  in  conjunction  with  the 
variational  multiscale  formulation  it  outperforms  the  current  state  of  art. 

The  results  of  this  research  have  been  published  in  refereed  journals  in  the  following 
articles  [16-19]. 


2 


Chapter  2 

Variational  Germano  Identity 


2.1  Introduction 

Consider  the  weak  form  of  an  abstract  (possibly  non-linear)  partial  differential  equation, 
viz.,  find  u  G  V,  such  that 

B(w,u)  =  (w,  f),  Vw  e  V.  (2.1) 

Here  B(-,  •)  :  V  x  V  — >  R  is  a  semi-linear  form  that  is  linear  in  its  first  slot,  V  is  the  space 
of  weighting  functions  and  trial  solutions,  and  /  €  L2(Q.)  is  the  prescribed  forcing  function. 
The  Galerkin  approximation  to  the  weak  form  is  given  by:  find  uh  €  Vh,  such  that 

B(wh,  uh )  -  (w\  /),  Vwh  e  Vh,  (2.2) 

where  Vh  C  V  is  a  finite  dimensional  subspace. 

In  several  applications  the  Galerkin  method  does  not  yield  good  numerical  solutions. 
In  such  cases,  other  numerical  methods  with  solutions  that  are  close  to  an  “optimal” 
representation  of  the  continuous  solution  in  Vh,  are  derived  (see  for  example  [20,21]).  That 
is  methods  with  solutions  uh, 


uh  «  vh  =  PV  (2.3) 

where  vh  is  the  optimal  representation  of  u  in  Vh,  and  :  V  — *  Vh  is  an  appropriate 
mapping.  Several  definitions  of  P*  and  hence  vh  may  be  considered.  For  example,  in  a  finite 
element  context,  vh  may  be  the  nodal  interpolant,  or  the  L2  or  H 1  projection  of  u  on  to  Vh. 
A  large  class  of  such  numerical  methods  may  be  formally  expressed  as:  find  uh  £  Vh,  such 
that 


B(wh,uh)  +  M{wh,uh,f;h,c)  =  (wh,f),  Vwh  6  Vh.  (2.4) 

Here  Mh( •,  •,•;•)•)  :  V  x  V  x  L2(Q)  x  R  x  Ip  ->  E  is  the  model  term,  which  is  a  functional  of 
the  weighting  function  wh,  the  trial  solution  uh,  and  the  forcing  function  /.  It  also  depends 
on  the  grid  or  mesh  size  h,  and  a  vector  of  parameters  c  =  [ci,  ■  ■  ■  ,  cp]  6  Rp.  Note  that  the 
solution  of  (2.4)  is  different  form  the  Galerkin  solution  (solution  of  (2.2)). 
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In  this  manuscript  we  consider  a  numerical  method  (such  as  the  one  in  (2.4))  with 
unknown  parameters  c,  and  derive  a  methodology  for  determining  these  dynamically.  This 
is  achieved  by  requiring  the  solution  of  the  numerical  method  to  be  equal  to  the  optimal 
representation  of  the  continuous  solution  on  Vh  and  its  subspaces.  The  result  is  an  expression 
for  c  in  terms  of  uh,  f,  the  mesh  size  associated  with  Vh  and  its  subspaces,  and  the  analog  of 
the  operator  Fh  on  each  subspace.  Similar  methods  for  computing  numerical  parameters  have 
been  developed  by  several  researchers  (see  for  example  [22-25]).  However,  these  methods 
typically  involve  either  analytical  or  numerical  solution  of  a  simplified  version  of  the  original 
or  an  auxiliary  (usually  the  adjoint)  PDE.  In  contrast  to  these,  the  proposed  method  does 
not  require  such  a  solution.  Instead  it  requires  restrictions  of  the  numerical  solution  on  to 
coarse  function  spaces.  These  restrictions  are  easily  computed  and  do  not  add  significantly 
to  the  overall  computational  cost  of  the  method. 

The  starting  point  of  our  development  is  the  variational  counterpart  of  the  Germano 
identity.  The  Germano  identity  is  a  popular  tool  for  computing  the  magnitude  of  the  eddy 
viscosity  in  the  large  eddy  simulation  of  turbulent  flows.  It  was  initially  derived  for  the 
filtered  Navier-Stokes  equations  in  [15].  In  this  manuscript  we  propose  a  generalization  of  the 
variational  form  of  this  identity,  and  describe  how  it  may  be  used  to  design  better  numerical 
methods.  We  first  apply  it  to  the  linear  advection-diffusion  equation  to  dynamically 
evaluate  the  mesh-dependent  diffusivity  for  the  finite  element  approximation  of  this  equation. 
Thereafter  we  apply  it  to  evaluate  the  Smagorinsky  coefficient  in  the  large  eddy  simulation 
(LES)  of  the  decay  of  homogeneous  isotropic  turbulence.  In  both  cases  the  resulting 
numerical  method  is  found  to  perform  well.  For  the  advection-diffusion  equation,  the  results 
compare  favorably  with  the  Streamline  Upwind  Petrov-Galerkin  (SUPG)  method,  which 
produces  nodally  exact  solutions.  For  the  turbulent  flow  problem,  the  results  are  generally 
more  accurate  than  the  constant-coefficient  and  the  traditional  dynamic  Smagorinsky 
models.  It  is  remarkable  that  the  same  formulation  leads  to  accurate  methods  for  these 
two  significantly  different  problems. 

The  format  of  this  chapter  is  as  follows:  In  Section  2,  we  describe  the  variational 
Germano  identity  and  derive  its  extension.  In  Section  3,  we  use  this  extension  to  derive 
an  expression  for  the  unknown  parameters  in  a  numerical  method,  which  is  the  main 
result  of  this  chapter.  In  Section  4,  we  apply  this  result  to  the  linear  advection  diffusion 
to  derive  a  non-linear  dynamic  diffusivity  method ,  and  compare  the  performance  of  this 
method  with  the  streamline-upwind  Petrov-Galerkin  (SUPG)  method.  In  Section  5  we 
apply  the  same  methodology  to  compute  the  Smagorinsky  coefficient  in  the  LES  of  the 
decay  of  homogeneous  isotropic  turbulence  in  three  dimensions,  and  compare  our  results 
with  the  constant-coefficient  and  the  traditional  dynamic  Smagorinsky  models.  We  end  with 
concluding  remarks  in  Section  6. 


2.2  The  Variational  Germano  Identity 

In  the  context  of  the  filtered  Navier-Stokes  equations  the  Germano  identity  is  derived  in  [15]. 
In  this  section  we  apply  its  variational  counterpart  to  the  abstract  PDE  given  by  (2.4). 

We  begin  by  requiring  the  solution  of  the  numerical  method  ( uh )  to  be  equal  to  vh. 
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Recall  that  vh  =  Fhu,  is  the  optimal  representation  of  Win  Vh.  Setting  uh  =  vh  in  (2.4), 

M(wh,vhJ-,h,c)  =  -(B(wh,vh)-(whJ)y  VwheVh.  (2.5) 

This  equation  may  be  considered  as  a  condition  on  the  model  term  that  is  necessary  for  the 
numerical  solution  to  be  equal  to  the  optimal  representation  of  the  continuous  solution  (that 
is  uh  =  vh).  Now  consider  the  finite  dimensional  subspace  V^1  C  Vh.  Let  vhl  —  Fhlu  be 
the  optimal  representation  of  u  in  Vhl  and  Fhl  :  V  — ►  V^1  be  the  appropriate  mapping.  We 
consider  the  following  numerical  method  which  is  obtained  by  replacing  h  with  hi  in  (2.4). 
Find  uhl  G  V'11 ,  such  that 

B(wh\uhl)  +  M(whltukltf]hi,c)  =  (wh\f),  Vwhl  €  Vhl.  (2.6) 

We  assume  that  the  solution  of  this  numerical  method  is  equal  to  the  optimal  representation 
of  it  in  Vhl .  That  is,  the  same  functional  form  of  the  model  that  leads  to  the  optimal  solution 
on  Vh,  also  leads  to  the  optimal  solution  on  Vhl.  Thus,  requiring  uhl  =  vhl  in  (2.6) 

M (whl ,vhl,  f]hi,c)  =  —  (jB(whl ,vhl)  —  (whl ,  /) Y  Mwhl  e  Vhl.  (2.7) 

Subtracting  (2.5)  from  (2.7)  we  arrive  at  the  variational  counterpart  of  the  Germano  identity, 
viz. 

M{wh\vk\f-huc)-M(wh\vh,f-h,c )  =  -(B(wh\vh')-B(whl,vh)y 

Vwhl  e  Vhl.  (2.8) 

Note  that  (2.8)  holds  for  weighting  functions  chosen  from  the  intersection  of  the  spaces  of 
weighting  functions  for  (2.5)  and  (2.7),  that  is  Vh  f)  V^1  =  Vhl.  Equation  (2.8),  which  may 
be  interpreted  as  a  consistency  condition  on  model  terms  at  two  different  scales,  involves 
only  the  finite  dimensional  representations  vh  and  vhl  of  the  exact  solution  u  and  not  u 
itself. 

Remark.  The  model  term  may  be  viewed  as  a  sub-grid  model  by  replacing  the  second  term 
on  the  right  hand  side  of  (2.5)  with  B(wh,  u).  This  yields 

M(wh,vh,f]h,c)  =  -(^B{wh,vh)-B{wh,u)j,  VwheVh.  (2.9) 

The  right-hand  side  of  (2.9)  represents  the  effect  of  that  part  of  u  which  is  not  contained  in 


2.2.1  Extension  of  the  Variational  Germano  Identity 

In  some  circumstances  (for  models  with  more  than  one  parameter)  the  variational  Germano 
identity  may  not  yield  sufficient  relations  to  determine  the  parameters  of  a  model.  For  these 
cases  we  need  to  generalize  (2.8). 

Consider  a  hierarchy  of  finite  dimensional  function  spaces  Vhj  C  Vhj~l  C  •  •  •  C  V^2  C 
Vhl  C  Vh.  For  a  given  space  Vh] ,  j  G  N(l,  J),  let  vhj  =  Fh^u  be  the  optimal  representation 
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of  the  continuous  solution  and  Fhj  :  V  — »  Vhj  be  the  appropriate  mapping.  Following  the 
procedure  described  in  the  previous  section,  and  assuming  that  the  same  functional  form  of 
the  model  yields  optimal  solutions  for  the  numerical  method  posed  on  all  subspaces  of  Vh, 
we  arrive  at  the  generic  Germano  identity  viz, 

M (whj ,  vhj ,  /;  hj,c)  —  M(whj  ,vh,f]h,c)  =  -  (B(wh> ,  vhj )  -  B(whj ,  vh)  ) , 

Vwhj  e  vh\  j  =  i,--- ,42.10) 

Note  that  for  a  given  j,  (2.10)  holds  for  weighting  functions  chosen  from  the  intersection  of 
the  spaces  Vh  f)  Vhj  =  Vhj . 

2.3  Application  Of  The  Variational  Germano  Identity 

In  this  section  we  describe  how  (2.10)  may  be  used  to  derive  an  expression  for  c. 

For  arbitrary  weighting  functions  the  variational  Germano  identity  (2.10)  represents  a 
relation  that  involves  the  functions  vh,f,  vhl ,  •  •  •  ,  vhj ,  and  the  parameters  h  and  c.  Here 

h  =  (2-11) 

The  dependence  on  vhl,  •  ■  •  ,  vhj  may  be  eliminated  by  expressing  vhj  in  terms  of  vh.  This 
is  possible  if 

pfyp*  =  pfy,  (2-12) 

since  this  implies  vhj  =  Fhjvh.  Assuming  Vhj  C  Vh,  it  can  be  verified  that  (2.12)  is  satisfied 
at  least  for  the  following  cases: 

1.  When  Ph  and  P,ij  are  interpolation  operators. 

2.  When  Vh  C  Hm(Cl),  and  and  are  //’"(fl)-projections  from  V  to  Vh  and  Vhj 
respectively,  with  n  <m. 

Using  (2.12),  (2.10)  is  reduced  to 

M(whj ,  Fhjvh,  /;  hj,  c )  -  M(w h* ,  vh,  /;  h,c)  =  -  (B{wh^ ,  P^  v'1)  -  B(w^ ,  V)) , 

Vwh>  €  Vh\  j  =  l,---  ,J.(2.13) 

In  this  form,  for  an  arbitrary  weighting  function,  the  variational  Germano  identity  represents 
a  relation  that  involves  the  functions  vh  and  /,  the  parameters  h  and  c,  and  the  set  of 
operators  P,  given  by 

P  =  {P'V--  ,P/lJ}-  (2-14) 

Let  Nhj  =  dim(V^).  Then  for  every  j  (2.13)  represents  Nhj  relations.  Experience  with 
the  Germano  identity  has  shown  that  these  relations  should  be  interpreted  in  a  global  sense. 
This  may  be  accomplished  by  either  of  the  following  methods. 
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Dissipation  method  This  approach  involves  choosing  wh]  =  vhj  =  Ph*vh  in  (2.13).  This 
yields 

M(Ph>vh,  Fhjvh,  /;  hj,c)  -  M(Phjvh,vh,  /;  h,c)  =  -(^B(Fhjvh,  Ph>vh)  -  B(Fhjvh,  vh)j, 

j  =  (2-15) 

Equation  (2.15)  is  motivated  by  previous  work  in  the  context  of  filtered  Navier-Stokes 
equations  [15].  In  that  case  the  model  contribution  in  (2.15)  (that  is  the  left  hand  side 
of  the  equation)  represents  the  difference  in  the  dissipation  of  the  total  turbulent  kinetic 
energy  induced  by  LES  models  acting  at  two  different  scales.  For  this  reason  we  refer  to  this 
approach  as  the  dissipation  method. 

Equation  (2.15)  represents  a  set  of  J  scalar  equations  for  the  P  parameters  c  = 
[ci,  •  •  •  ,  cp]  that  involve  vh,  /,  h  and  P.  When  J  =  P,  these  equations  may  be  solved 
to  determine  the  model  parameters.  Assuming  this  is  feasible  the  result  is  a  relation  of  the 
form 


c  =  n(vh,f;h]P).  (2.16) 

Least-squares  method  An  alternate  approach,  which  is  motivated  by  the  work  of  Ghosal 
et  al.  [26]  and  Lilly  [27]  involves  selecting  wkl  as 

wh*  =  <t>hj{x)  (2.17) 

where  4*  a  ,  A  =  1,  •  •  •  ,  Nhj  are  functions  that  span  Vhj .  Using  (2.17)  in  (2.13), 

M{4>hJ(,Fh^vh,  f  ; hj, c)  -  M(4>hj,vh, /; h, c)  =  -(B(<f>%,Ph*vh)  -  B(cfy ,‘u'*)), 

A=1,..-,W\  j  =  l,---,J. 

(2.18) 

The  equation  above  represents  N  —  J2j=i  Nhj ,  scalar  relations.  An  expression  that  is 
formally  identical  to  (2.16)  may  be  obtained  from  this  equation  by  finding  the  parameters  c 
that  minimize  the  square  of  the  residual  of  these  relations. 

Once  an  expression  for  7r  is  derived  using  either  the  dissipation  or  the  least-squares 
method,  the  numerical  method  may  be  written  as 

B{wh,uh)  +  M(wh,uh,f-,h,c )  =  ( wh,f ),  Vwh  e  Vh  (2.19) 

c  =  7r(u\/;/i;P),  (2.20) 

By  construction,  this  method  has  one  of  the  following  special  properties: 

1.  When  the  parameters  are  determined  using  the  dissipation  method,  a  subset  of 
the  conditions  that  are  necessary  for  the  solution  of  the  numerical  method  on 
Vh,  Vhl,  •  •  •  ,  Vhj  to  be  equal  to  the  optimal  representation  of  the  continuous  solution, 
are  satisfied. 
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2.  When  the  parameters  are  determined  using  the  least-squares  method,  the  conditions 
that  are  necessary  for  the  solution  of  the  numerical  method  on  Vk,  V^1,  •  •  •  ,  Vhj  to 
be  equal  to  the  optimal  representation  of  the  continuous  solution,  are  satisfied  in  a 
least-squares  sense. 

The  following  additional  remarks  may  be  made  about  the  numerical  method  (2.19)-(2.20). 
Remarks 

1.  During  a  numerical  simulation,  each  quantity  that  appears  in  (2.20)  is  known,  hence  this 
expression  closes  the  numerical  method. 

2.  7r  is  typically  non-linear  in  uh,  and  hence  the  resulting  numerical  method  is  also  non¬ 
linear. 

3.  The  fact  that  7r  depends  on  P,  which  in  turn  depends  on  the  definition  of  the  optimal 
solution,  indicates  that  different  numerical  methods  are  obtained  for  different  choices  of  the 
optimal  solution. 

4.  Evaluating  7r  using  either  the  dissipation  or  the  least-squares  method,  and  hence 
implementing  the  proposed  numerical  method,  does  not  require  the  analytical  solution  of 
the  continuous  problem.  Thus  this  method  may  be  used  for  solving  complex  non-linear 
partial  differential  equations  where  such  solutions  are  unavailable. 

5.  In  the  development  above,  for  simplicity,  we  have  assumed  that  c  is  constant  in  fh  When 
this  is  not  the  case  a  piecewise-constant  approximation  to  c  may  be  used  and  the  procedure 
described  above  repeated  on  individual  subdomains  where  c  is  constant. 


2.4  Linear  Advection  Diffusion  Equation 

In  this  section  we  apply  the  methodology  developed  in  the  previous  section  for  solving  the 
linear  advection-diffusion  problem.  For  this  equation  a  model  term  that  yields  nodally  exact 
solutions  in  one  dimension  can  be  derived  analytically.  Hence  the  results  of  (2.19)-(2.20) 
may  be  compared  with  this  “ideal”  model. 

2.4.1  Problem  Description 

We  consider  the  advection  diffusion  equation  in  f2  =  (0, 1),  with  homogeneous  Dirichlet 
boundary  conditions  and  forcing  /.  The  strong  form  of  the  problem  is  given  by 

Cu  =  -uutXX  +  auiX  =  /,  x  e  (0, 1)  (2-21) 

u(0)  =  u(l)  =  0.  (2.22) 

where  v  is  the  diffusivity,  assumed  constant  in  D,  a  =  1,  is  the  prescribed  velocity,  and  /  is 
the  forcing  function  given  by 

'(*)  =  <*  -*»>  =  {*-  ItZ  (2'23) 

We  have  experimented  with  several  values  of  xo  and  have  obtained  results  that  are  insensitive 
to  it.  As  a  representative  value  we  choose  x0  =  1/2. 
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An  equivalent  weak  or  variational  form  of  (2.21)  and  (2.22)  is  given  by  (2.1),  where 


B(w,  u ) 

K/) 


v{w,x,u,x)  -  a(wiX,u), 


and  V  =  {v\v  E  if^Q),^)  =  v(l)  =  0}. 


(2.24) 

(2.25) 


2.4.2  Variational  Germano  Identity 

Preliminaries 

Model  term  The  finite  dimensional  approximation  of  (2.21)-(2.22)  inclusive  of  a  model 
is  given  by  (2.4).  For  the  model  we  choose  a  residual-based  term  given  by 

M(w\  u\  /;  h,  c )  =  Clhc\w%,  Cuh  -  f){ j,  (2.26) 

where  12  is  the  union  of  all  element  interiors,  h  is  the  element  size  of  the  finite  element 
discretization,  and  c  =  {01,02}  are  the  parameters  that  will  be  determined  using  (2.15). 
Since  the  number  of  unknown  parameters  P  =  2,  we  select  the  number  of  subspaces  of  V7*, 
denoted  by  J  =  2  also. 


Function  spaces  To  solve  (2.4)  we  use  a  uniform  mesh  of  linear  finite  elements.  With 
this  specification  the  model  term  (2.26)  reduces  to 

M(w\  u\  /;  h,  c)  =  Clhc>(w%  uhx  -  /).  (2.27) 

The  spaces  Vh,  V711  and  Vh2  are  constructed  from  shape  functions  associated  with  a  uniform 
mesh  of  linear  finite  elements  of  size  h,  2h  and  4 h  respectively.  Thus  h  =  {h,2h,4h}.  It 
easily  verified  that  this  choice  satisfies  Vhi  C  Vhl  CVh.  ■ 

Operators  Since  J  —  2,  the  set  P  =  (P7*1,  P712}.  The  operators  P/l,  P7*1  and  P7*2  are  chosen 
to  be  the  nodal  interpolation  operators  at  scales  h,  2 h  and  4 h  respectively.  Note  that  this 
choice  satisfies  (2.12). 

Determining  7r 

We  now  specialize  (2.15)  to  the  advection-diffusion  problem.  Using  (2.27)  in  (2.15), 
recognizing  that  £?(•,  •)  is  a  bilinear  form,  and  that  for  our  choices  of  Vhj  and  P,lj , 

(^,«t-/)  =  (^,(Pft74-/).  V/eV^,  J  =  1,2,  (2.28) 

we  arrive  at 

ClhC2((hj/hr  -  1)  =  F(vh,  j  =  1,2,  (2.29) 
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where 


F(vh  B(Fhwh,vh  -Fh^vh) 

K  }-((FhJv^x,(F  WO,*-/)' 

(2.30) 

Solving  (2.29)  for  c\  and  c2  yields 

C2  -  7T2(^h,  /;  /j.;  P)|/,.={^2^4M  —  7T2(^,  /;  /i;  P)  =  In 

l)(ln2)-(2.31) 

Cl  -  tfiOA  /;  h ;  P)\h={Wh}  ~  /;  h;  P)  -  ^ . 

(2.32) 

2.4.3  Numerical  Solution 

Nonlinearity  The  numerical  method  (2.19)-(2.20)  specializes  to 

B(wh,  uh )  +  vh{whx,  uhx  -f)  =  ( wh ,  /),  Vwh  €  Vh 

(2.33) 

»h  =  g\u\f-h-  P) 

=  al(uhJ;h]P)h^uh'f'h’^ 

(2.34) 

where  £?(•,  •)  is  defined  in  (2.24),  and  7^  and  7r2  are  given  by  (2.32)  and  (2.31).  We  term 
this  method  the  dynamic  diffusivity  method ,  since  the  numerical  diffusivity  is  not  determined 
a-priori,  but  is  determined  based  on  the  solution. 

To  solve  this  nonlinear  problem  we  introduce  an  artificial  “time”  variable  t,  and  instead 
of  (2.34)  we  solve 

di/h 

~  =  gh{u\f-h-  P)-uh,  (2.35) 

with  vh(Q)  =  0.  The  steady  state  solution  of  this  equation  is  the  solution  to  (2.34).  We  use 
the  following  simple  scheme  to  discretize  this  equation  in  time. 

A*  +  A()  =  YTKt^(t)  +  TTSls‘(“',(t)’ h ' p)l  (2'36) 

where  uh(t )  is  solution  to  (2.33).  This  scheme  corresponds  to  a  forward  Euler  step  in  the 
non-linear  term  and  a  backward  Euler  step  in  the  linear  term.  We  have  used  the  absolute 
value  of  gh  in  this  equation,  in  order  to  avoid  the  scheme  from  diverging  whenever  (rarely) 
a  negative  value  of  gh  is  encountered.  Similar  precautions  are  taking  when  evaluation  eddy 
viscosities  in  turbulence  models  [15,26].  For  the  examples  shown  below  we  choose  A t  =  3/7. 

Results  The  diffusivity  that  renders  the  numerical  solution  nodally  exact  is  given  by 

=  yC.  (2-37) 

where  the  non-dimensional  parameter  £e  is  given  by 

f  =  coth  a--,  (2.38) 

a 
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and  a  =  is  the  mesh  Peclet  number.  Using  this  diffusivity  in  (2.33)  leads  to  the  SUPG 
method  [28]. 

First,  we  study  the  convergence  of  the  dynamic  diffusivity  uh,  obtained  by  solving  (2.33) 
and  (2.35)  to  ve.  In  particular  we  compare  with  £e,  where  £h  is  given  by 


„  sen 

ah/2  ’ 


(2.39) 


where  uh(T)  is  the  converged  or  the  steady-state  value  of  uh.  The  results  of  this  comparison 
are  shown  in  Figure  2.1  for  a  6  (3.2  x  10-2, 3.2  x  103).  For  all  results  h  =  1/52.  We  observe 
that  the  values  of  t/1  and  £e  are  in  agreement  in  both  the  advective  and  diffusive  limits 
(a  — >  oo  and  a  — >  0  respectively).  There  are  some  noticeable  differences  in  the  diffusive 
limit  with  dynamic  diffusivity  being  slightly  larger  than  the  analytical  value.  However,  it  is 
worth  noting  that  in  this  limit  the  numerical  diffusivity  is  much  smaller  than  the  molecular 
diffusivity  (same  for  every  method),  which  dominates  the  solution. 


Figure  2.1:  Variation  of  the  non-dimensional  viscosity  for  SUPG  (£e)  and  the  dynamic 
diffusivity  method  (/h)  as  a  function  of  the  mesh  Peclet  number  (a). 

In  Figures  2.2  through  2.5,  we  have  plotted  the  Galerkin,  the  SUPG  (which  is  also  the 
nodal  interpolant  of  the  exact  solution)  and  the  dynamic  diffusivity  solutions  for  different 
values  of  mesh  Peclet  number.  Figure  2.2  corresponds  to  a  value  of  a  =  3.2  x  10”2.  We 
observe  that  all  solutions  including  the  Galerkin  solution  are  very  accurate.  In  Figure  2.3  we 
have  plotted  the  solutions  corresponding  toa=  1.0.  We  can  now  observe  differences  in  the 
Galerkin  solution  and  other  solutions.  There  is  a  bump  in  the  Galerkin  solution  at  x  =  0.97 
which  is  absent  from  the  SUPG  (the  nodal  interpolant)  and  the  dynamic  diffusivity  solution. 
In  Figure  2.4  we  have  plotted  the  solutions  for  a  —  3.2  x  103.  This  represents  the  advective 
limit.  We  observe  that  the  Galerkin  solution  has  large  spurious  oscillations  that  are  absent 
from  the  SUPG  and  dynamic  diffusivity  method.  In  order  to  compare  the  dynamic  method 
with  the  SUPG  method,  we  have  plotted  these  solutions  in  Figure  2.5  without  the  Galerkin 
solution.  We  observe  that  the  solutions  are  virtually  indistinguishable. 
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Figure  2.2:  Galerkin,  SUPG  and  dynamic  diffusivity  solutions  for  a  =  3.2  x  10  2. 

In  Figure  2.6,  we  present  the  convergence  of  the  diffusivity  of  the  dynamic  method  as  a 
function  of  iteration  number,  for  different  values  of  Peclet  numbers.  For  all  cases  we  start 
with  no  diffusivity  (that  is  the  Galerkin  method)  and  within  about  10  iterations  converge  to 
the  SUPG  method. 


2.5  Incompressible  Navier-Stokes  Equation 

In  this  section  we  apply  the  variational  Germano  identity  to  the  more  challenging  problem 
of  solving  turbulent  flows  governed  by  the  incompressible  Navier-Stokes  equations  on  grids 
for  which  the  mesh  size  is  much  larger  than  the  dissipation  length  scale. 


2.5.1  Problem  Description 

We  consider  the  problem  in  Q  =  flx]Ti,  7a[,  where  Q  =  [0, 27r]3  is  the  spatial  domain  and 
]Ti,  T%[  is  the  time  period  of  interest.  The  strong  form  of  the  problem  is  given  by 


f  u<t  +  V  •  (u  <g>  u)  +  Vp  -  z/V2«  =  0,  in  Q 
\  V  •  u  =  0,  in  Q 


(2.40) 


In  the  above  equations  U  =  [u,p]T,  where  u  is  the  fluid  velocity  field  and  p  is  the  pressure. 
In  addition,  u  is  the  kinematic  viscosity  and  the  symbol  <g>  denotes  the  outer  product  of  two 
vectors.  The  boundary  of  Q,  is  denoted  by  <9f 2.  It  is  comprised  of  six  faces,  Fj(0),  Fj(27r), 
j  —  1, 2, 3,  where 


r,-(c)  =  {sc  €  dQ  |  Xj  =  c} 
It  is  assumed  the  solution  is  periodic,  that  is 


(2.41) 


u(x  +  27rej,£)  =  u(x,t),  x  G  Fj(0),  t  gJTx, 


(2.42) 
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Figure  2.3:  Galerkin,  SUPG  and  dynamic  diffusivity  solutions  for  a  —  1. 

where  ej  is  the  Cartesian  basis  vector  in  the  Xj  direction.  The  initial,  divergence-free  velocity 
field  is  given  by 

u(x,  0)  =  u0(x).  (2.43) 

When  seeking  approximate  numerical  solutions  to  the  Navier-Stokes  equations  one  may 
either  (a)  develop  and  approximate  a  spatial  variational  formulation  of  the  problem  to 
arrive  at  a  set  of  coupled  ordinary  differential  equations  and  solve  these  using  a  standard 
time  marching  scheme,  or  (b)  discretize  the  space- time  domain  into  slabs,  and  develop  and 
approximate  a  space-time  variational  formulation  using  the  discontinuous  Galerkin  method, 
leading  to  a  set  of  coupled  algebraic  equations.  Methods  based  on  approach  (a)  are  commonly 
referred  to  as  semi-discrete  methods,  and  those  based  on  approach  (b)  are  referred  to  as  fully 
discrete  space-time  methods.  In  our  work  we  consider  the  more  commonly  used  semi-discrete 
method.  A  direct  consequence  of  this  choice  is  that  numerical  method  developed  in  Section  2 
(which  is  directly  applicable  to  space-time  methods)  must  be  modified.  These  modifications 
are  described  in  the  following  paragraphs. 

The  weak  formulation  (in  tt)  of  (4.1),  (4.30)  and  (4.3)  is  given  by  (see  [29]  for  example): 
Find  u(-,  t)  €  V,  such  that 

B{w,u)  =  0,  VwG  V,Vt  GjTi.Tat,  (2.44) 

where  the  semi-linear  form  is  defined  as 

B(w,u )  =  (w,  uit)  —  (Vto,  u  <g>  u)  -I-  2v(Vsw,  \7su)  +  (V  •  m,  V~2(V  •  V  •  u  ®  t(]B)45) 

Note  that  the  pressure  and  the  divergence  terms  are  absent  from  the  definition  of  the  semi- 
linear  form.  This  is  because  the  divergence-free  constraint  for  the  velocity  field  is  built  into 
the  definition  of  the  function  space  V, 

V  =  jv|o  =  VkeikX’  v*-k  =  +  lfc!2)  <  °°}-  (2-46) 

ke  7?  kez3 
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Figure  2.4:  Galerkin,  SUPG  and  dynamic  diffusivity  solutions  for  a  =  3.2  x  103. 

In  the  equation  above,  v ^  is  the  Fourier-coefficient  of  v,  k  is  the  corresponding  wave- vector, 
and  the  superscript  *  denotes  the  complex-conjugate  of  a  quantity.  The  terms  within  the 
curly  brackets  imply  that  vector  fields  contained  in  V  are  real- valued,  divergence-free  and 
belong  to  H1  (Q)  (the  energy  and  the  enstrophy  of  the  flow  are  finite  for  t  e]Tj,  T2[).  The 
periodic  boundary  conditions  are  also  built  into  the  definition  of  V. 

2.5.2  Variational  Germano  Identity 

Model  term  The  finite  dimensional  approximation  of  (4.11)  and  (4.3),  inclusive  of  a 
model,  is  given  by:  Find  uh(-,t )  E  Vh  such  that 

B(wh, uh)  +  M(wh,  uh\  h, c)  —  0,  Vwh  eVh,VtE}TuT2[.  (2.47) 

For  the  model  we  choose  the  Smagorinsky  eddy  viscosity  model  [30]  given  by, 

M(wh ,  uh ;  h,  c)  =  2(clh)2(Wsw\  | Vsuh\  VV1),  (2.48) 

where  h  is  the  grid  (or  mesh)  size,  and  Ci(t)  is  the  time- dependent  viscosity  coefficient  that 
will  be  determined  using  the  methodology  developed  in  Section  2. 

For  the  semi-discrete  approximation  of  the  Navier-Stokes  equations  considered  herein, 
the  variational  Germano  identity  may  be  derived  using  arguments  that  are  analogous  to 
those  used  in  Sections  2  and  3.  The  end  result  is  (2.13),  which  now  holds  Vt  g]Ti, T2[.  The 
application  of  the  dissipation  and  the  least  squares  method  to  this  equation  results  in  (2.15) 
and  (2.18)  respectively,  which  now  hold  Vt  E]Xj,  T2[.  Consequently,  these  equations  may 
be  used  to  determine  the  time-dependent  parameters  in  the  model  term.  In  the  following 
development,  for  the  Smagorinsky  model,  we  use  (2.15)  to  determine  c\  as  a  function  of  time. 

Function  spaces  Since  the  number  of  unknown  parameters  is  one  (P  =  1),  we  select  the 
number  of  subspaces  of  Vh  also  to  be  one  (J  =  1).  We  use  a  finite-dimensional  Fourier- 
spectral  basis  to  approximate  our  solution  in  O.  Formally,  the  spaces  Vh  and  V^1  are  given 
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Figure  2.5:  SUPG  and  dynamic  diffusivity  solutions  for  a  —  3.2  x  103. 


by 


V"  =  {u  EV;vje  =  0,  Ifcloo  >  kh}  (2.49) 

Vhl  =  {veV;vk  =  0,\  k\00>kh'},  (2.50) 

where  Ifcloo  =  max{|fci|,  |fc2|,  |^3|},  and  kh  =  ir/h,  and  khl  =  Tt/hy  are  the  cut-off 
wavenumbers.  In  particular,  we  choose  kh  =  16,  and  khl  =  8. 

Operators  Since  J  =  1,  the  set  P  =  The  operators  P/(  and  P/tJ  are  chosen  to 

be  the  Hl(¥l)  projections  of  functions  in  V  into  and  Vhl  respectively.  Note  that  these 
projections  satisfy  the  following  property. 

Let  w(-,  f)  G  V  be  represented  in  terms  of  the  infinite  series 

u(x,t)  =  ^Wfc(t)eifcx ,  (2.51) 

fc 

where  uk(t)  are  the  Fourier-coefficients  of  u.  Then  Fhu  is  given  by 

P  hu(x,t)  =  (2.52) 

|fc|oo  <kh 

and  likewise  P^w  is  given  by 

¥hlu{x,t)  =  J2  ^fcW6^'*-  (2-53) 

|fc|oo<fchl 

In  other  words,  these  projections  correspond  to  the  sharp-cutoff  filters  in  the  wavenumber 
space. 


15 


Figure  2.6:  Variation  of  uh  for  the  dynamic  diffusivity  method  as  a  function  of  iteration 
number. 


Determining  7r 

We  now  specialize  (2.15)  to  the  incompressible  Navier-Stokes  equation.  Using  the  definition 
of  the  model  (2.48)  and  the  semi-linear  form  (2.45)  in  (2.15),  we  have 

Ci  =  _ 

/l  {yvh\vh'  ®vh')  -  (Vvh',vh®vh)  ,  fn 

y  2hf(ysvh\\Vsvhi\Vsvhi)  -  h2{Vsvhi,\Vsvh\Vsvh)’  '  U  {  '  ’ 

where  vhl  =  Fhlvh.  In  deriving  (2.54)  we  have  utilized  the  fact  that  the  contribution  from 
terms  in  B(w,  u )  that  are  linear  in  u  is  zero  due  to  the  specific  form  of  the  basis  functions 
and  projection  operators  IP'1  and  Fhl .  We  have  also  made  use  of  the  fact  that  vh  and  vhl  are 
divergence  free. 


2.5.3  Numerical  Solution 

The  proposed  numerical  method  for  solving  the  incompressible  Navier-Stokes  equations  is 
given  by  (2.47)  and  (2.48).  The  Smagorinsky  parameter  that  appears  in  the  model  is  given 
by  (2.54)  with  vh  and  vhl  replaced  by  uh  and  uhl  respectively.  That  is 


l _ ( <8>  )  -  (' Vuh'  ,Uh®Uh) _  v  ]T  rp  r  /2«\ 

Cl  2h21(Vsuh',\Vsuhl\Vsuh')-h2{Vsuhl,\X7suh\X7sUh),  1  U  21,1  j 

where  uhl  =  Fhluh.  We  refer  to  the  method  implied  by  (2.47),  (2.48)  and  (2.55)  as  the 
variational  dynamic  Smagorinsky  method.  These  equations  represent  a  closed  system  of 
ODEs  which  may  be  integrated  in  time  to  yield  the  velocity  field  in  ]Ti,T2[. 
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Results  To  test  the  performance  of  the  proposed  numerical  method  we  generate  a  well- 
resolved  benchmark  solution  on  a  2563  mesh.  For  this  solution  we  use  a  random  initial 
condition  with  an  energy  spectrum  (energy  density  in  wavenumber  space)  given  by  E(k )  = 
(q2 /2A)(k4/kp)  exp(— 2(k/kp)2),  where  the  initial  turbulent  kinetic  energy,  ^-  =  |,  and  the 
spectrum  is  peaked  at  kp  =  3.  The  same  spectrum  was  used  in  [31]  to  study  the  decay  of  low 
Reynolds  number  homogeneous  isotropic  turbulence.  The  phase  of  the  modes  for  the  initial 
velocity  field  is  chosen  randomly  and  each  mode  satisfies  the  divergence-free  condition.  The 
solution  is  allowed  to  evolve  according  to  the  Navier-Stokes  equations  till  a  spectra  with  a 
physical  k~5/s  range  is  obtained  (see  Figure  2.7).  This  corresponds  to  t  —  Ti  «  2.2.  At 
this  time  the  Taylor-microscale  Reynolds  number  is  ~  90.  This  state  of  the  flow  is  chosen 


Figure  2.7:  Energy  spectra  E(k),  for  the  DNS  at  initial  and  final  time  (Ti&T2). 

as  an  initial  condition  for  all  numerical  solutions.  In  addition  to  the  variational  dynamic 
Smagorinsky  method  on  a  323  grid,  the  following  numerical  methods  are  included  in  the 
comparison. 

1.  A  well  resolved  Direct  Numerical  Simulation  (DNS)  solution  on  a  2563  mesh.  This 
corresponds  to  the  Fourier-Galerkin  method  applied  to  the  Navier-Stokes  equations 
(that  is  (2.47)  with  the  model  term  set  to  zero).  The  resolution  is  chosen  such  that  the 
largest  wavenumber  is  sufficient  to  capture  the  dissipation  or  the  Kolmogorov  length 
scale.  This  solution  is  treated  as  the  benchmark  solution.  Note  that  the  number 
of  degrees  of  freedom  in  the  well  resolved  DNS  is  512  times  more  than  the  other 
simulations. 

2.  An  under-resolved  (coarse)  DNS  solution  on  a  323  mesh.  This  method  is  also  given 
by  (2.47)  with  the  model  term  set  to  zero.  However  a  coarse  323  mesh,  on  which  the 
dissipation  length  scale  is  not  represented,  is  utilized.  This  solution  represents  the 
coarse  Fourier-Galerkin  approximation  of  the  Navier-Stokes  equations. 

3.  A  constant  coefficient  (static)  Smagorinsky  model  on  a  32s  mesh.  This  numerical 
method  is  given  by  (2.47)  and  (2.48),  where  the  Smagorinsky  constant  is  set  to 


17 


Ci  =  0.16,  a  value  obtained  by  assuming  an  infinite  inertial  range  (see  [32]).  The 
value  of  the  coefficient  is  not  changed  during  the  simulation. 

4.  The  traditional  dynamic  Smagorinsky  model  on  a  323  mesh.  This  numerical  method  is 
given  by  (2.47)  and  (2.48).  In  (2.48),  the  Smagorinsky  coefficient  is  calculated  using  the 
Germano  identity  [15]  applied  to  the  filtered  Navier-Stokes  equations.  The  coefficient 
is  given  by 


/! _ (Vnfe,  uh'  ®uh')~  (Vttfei ,  uh  ®  uh ) _  f  . 

1  y  2hl(Vsuh,  |Vsu/ll|Vsw/l1)  -  h?(Vsuh\  |V5^|VS^)’  J  15  2['  {  } 

Note  that  in  the  above  expression  for  ci,  uh  appears  in  the  weighting  function  slot 
of  the  second  term  in  the  numerator  and  the  denominator.  Whereas  in  (2.55),  which 
is  derived  from  a  variational  standpoint,  uhl  appears  in  these  slots.  In  the  following 
section  we  will  observe  that  this  difference  leads  to  a  smaller  value  for  ci  when  it  is 
calculated  using  (2.55). 

All  the  numerical  methods  described  above  lead  to  a  set  of  coupled  ODEs.  These 
equations  are  advanced  in  time  using  a  mixed  integration  scheme  in  which  the  exact 
integrating  factor  is  used  for  the  viscous  term  and  a  third-order  Runge-Kutta  method  is 
used  for  the  remaining  terms  (see  for  example  [33]). 

In  Figure  2.8,  we  have  plotted  the  Taylor  microscale  Reynolds  number  (Rex)  as  a  function 
of  time  for  the  benchmark  DNS  solution.  We  observe  that  in  the  time  period  of  interest  (that 
is  t  e]2.2,4.4[),  Rex  falls  from  90  to  about  62. 


Figure  2.8:  The  Taylor  microscale  Reynolds  number  (Re a)  for  the  DNS  as  a  function  of  time. 

In  Figure  2.9,  we  have  plotted  the  coefficient  ci  as  a  function  of  time  for  the  static,  the 
traditional  dynamic  and  the  variational  dynamic  Smagorinsky  formulations.  We  observe 
that  the  coefficient  for  the  static  Smagorinsky  formulation  is  the  largest.  For  the  traditional 
dynamic  method  the  coefficient  starts  out  high  but  settles  to  a  lower  value  of  0.14.  The 
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Figure  2.9:  Variation  of  the  Smagorinsky  coefficient  (ci)  as  a  function  of  time. 


coefficient  for  the  variational  dynamic  method  is  the  smallest,  and  settles  to  a  value  of  about 

0.12. 

In  Figure  2. 10,  we  have  plotted  the  energy  spectra  for  the  benchmark  DNS  and  the  various 
numerical  solutions  at  the  final  time  (T2).  From  this  figure  we  observe  that  the  coarse  DNS 
solution  is  completely  incorrect  indicating  that  a  model  term  is  necessary  to  compute  a 
reasonable  approximation  of  the  well-resolved  DNS  solution.  The  coarse  DNS  exhibits  a 
cusp  at  large  wavenumbers  caused  by  the  abrupt  truncation  of  the  energy  cascade  at  the 
cut-off  wavenumber.  Further,  through  the  coupling  induced  by  the  quadratic  term,  this 
accumulation  of  energy  at  high  wavenumbers  has  led  to  an  increased  transfer  of  energy  from 
the  low  wavenumbers,  which  has  in  turn  resulted  in  an  under-prediction  of  the  energy  at  these 
wavenumbers.  This  error  in  the  coarse  DNS  solution  is  rectified  by  the  Smagorinsky  model. 
However,  for  the  static  coefficient  case  we  observe  that  the  model  has  overcompensated  to 
some  extent:  the  large  dissipation  in  the  fine  scales  has  led  to  a  drop-off  in  the  spectrum  at 
the  cutoff  wavenumber,  which  has  in  turn  caused  a  net  decrease  in  the  transfer  of  energy  from 
the  coarse  scales,  leading  to  a  pile-up  of  energy  at  the  these  scales.  The  traditional  dynamic 
Smagorinsky  formulation,  which  has  a  smaller  viscosity,  is  seen  to  improve  these  results 
somewhat.  The  most  accurate  spectrum  is  achieved  by  the  variational  dynamic  formulation, 
which  has  the  smallest  Smagorinsky  coefficient. 

In  Figure  2.11,  we  have  plotted  the  resolved  turbulent  kinetic  energy  in  the  flow  as  a 
function  of  time.  We  observe  that  with  the  exception  of  the  coarse  DNS  solution,  all  solutions 
match  the  DNS  data  quite  well.  In  particular,  the  variational  dynamic  Smagorinsky  and  the 
constant  coefficient  Smagorinsky  methods  are  very  accurate  at  initial  and  intermediate  times 
respectively. 

Finally  in  Figure  2.12,  we  have  plotted  the  resolved  enstrophy  (L2  norm  of  vorticity) 
as  a  function  of  time.  We  observe  that  the  coarse  DNS  solution  severely  over-predicts 
this  quantity,  while  both  the  constant  coefficient  and  the  traditional  dynamic  Smagorinsky 
solutions  under-predict  it.  The  variational  dynamic  formulation  is  once  again  the  most 
accurate. 


Figure  2.10:  Energy  spectra  E(k),  for  the  benchmark  DNS  and  the  numerical  methods  at 
the  final  time  t  —  TV 

2.6  Conclusions 

Using  an  extension  of  the  variational  counterpart  of  the  Germano  identity,  we  have  derived  a 
methodology  for  evaluating  the  parameters  of  a  numerical  method.  The  parameters  obtained 
using  this  approach  satisfy  certain  conditions  necessary  for  the  numerical  method  to  yield 
optimal  solutions  over  a  set  of  finite  dimensional  subspaces.  We  have  applied  this  approach 
to  calculate  the  numerical  diffusivity  required  for  nodally  exact  solutions  to  the  advection 
diffusion  equation  in  one  dimension.  We  have  found  that  the  resulting  diffusivity,  which 
is  a  non-linear  function  of  the  numerical  solution,  converges  to  the  value  obtained  using 
the  SUPG  analysis.  In  addition  we  have  applied  the  same  approach  for  computing  the 
Smagorinsky  coefficient  in  the  decay  of  homogeneous  isotropic  turbulence.  In  this  case  we 
have  found  that  the  resulting  model  produces  results  which  are  generally  more  accurate  than 
the  constant-coefficient  and  the  traditional  dynamic  Smagorinsky  models. 

It  is  interesting  to  enumerate  the  differences  between  the  advection-diffusion  and  the 
turbulence  problems  that  were  solved  in  this  study:  (1)  the  former  was  a  linear  PDE  while 
the  latter  was  a  system  of  non-linear  PDEs,  (2)  the  former  was  posed  in  one  spatial  dimension 
while  the  latter  was  posed  in  three  dimensions,  (3)  the  latter  exhibited  chaotic  solutions  while 
the  former  did  not,  (4)  the  former  was  solved  using  the  finite  element  method,  while  the 
latter  was  solved  using  a  spectral  method.  Given  these  differences  it  is  remarkable  that  the 
same  methodology,  namely  the  variational  Germano  identity,  could  be  employed  to  design 
accurate  numerical  methods  for  these  two  problems.  In  the  future  we  propose  to  explore  the 
application  of  this  identity  to  other  physical  systems. 


Figure  2.11:  Total  resolved  kinetic  energy  ( q 2)  as  a  function  of  time. 


Figure  2.12:  Total  resolved  enstrophy  as  a  function  of  time 


Chapter  3 

Variational  Germano  Identity  Applied 
to  the  Spectral  Discretization  of  a 
Conservation  Law 


3.1  Introduction 

In  this  chapter  we  develop  a  numerical  method  for  the  spectral  approximation  of  non-linear 
conservation  laws.  These  laws  describe  a  broad  range  of  physical  phenomena  which  include 
the  dynamics  of  gasses,  the  flow  of  traffic  and  the  propagation  of  shallow  water  and  nonlinear 
acoustic  waves.  In  all  these  systems  we  are  interested  in  cases  when  the  physical  viscosity 
(or  diffusivity)  is  small  or  zero.  In  the  small  viscosity  case,  the  solution  to  such  systems 
develops  local  regions  of  large  spatial  and  temporal  gradients  called  shocks.  The  width  of  a 
shock  reduces  with  reducing  viscosity,  and  in  the  limit  of  zero  viscosity  the  solution  becomes 
discontinuous.  In  fact,  in  this  limit  in  order  to  ensure  unique  solutions,  the  conservation 
law  must  be  supplemented  with  an  entropy  production  inequality  and  conditions  that  relate 
jumps  in  conserved  quantities  across  the  shock  [34,35]. 

For  small  viscosities,  the  standard  Fourier-Galerkin  approximation  to  non-linear 
conservation  laws  becomes  unstable  if  the  shock  width  is  smaller  than  the  grid  size.  For 
a  large  class  of  problems  the  computational  cost  of  employing  a  grid  which  is  fine  enough  to 
resolve  a  shock  is  prohibitive  and  as  a  result  this  method  finds  limited  application.  Further,  in 
the  limit  of  zero  viscosity,  even  with  sufficient  grid  refinement,  the  Fourier-Galerkin  solution 
does  not  converge  to  the  unique  “physical”  solution  which  satisfies  the  entropy  production 
inequality.  To  overcome  these  difficulties  associated  with  the  Fourier-Galerkin  method, 
several  methods  have  been  proposed.  A  large  proportion  of  these  methods  involve  appending 
to  the  Fourier-Galerkin  formulation  a  numerical  viscosity  term  (see  [36]  for  example).  We 
choose  to  classify  different  numerical  viscosity  based  methods  on  the  basis  of  the  equations 
in  which  the  viscosity  appears.  To  accomplish  this,  we  introduce  the  concept  of  the  coarse 
and  the  fine  scale  equations  of  a  numerical  approximation  as  follows. 

In  a  Fourier-Galerkin  method,  the  residual  of  the  original  partial  differential  equation  is 
weighted  by  a  Fourier  mode,  integrated  over  the  domain,  and  the  result  is  set  to  zero.  This 
leads  to  a  finite  number  of  ordinary  differential  equations  (ODEs),  which  may  then  be  solved 
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to  determine  the  coefficients  in  the  Fourier  expansion  of  the  numerical  solution  (see  [37]  for 
example).  Note  that  in  a  Galerkin  method  the  same  set  of  modes  is  used  for  the  weighting 
functions  and  the  trial  solution.  Given  a  set  of  modes  that  comprises  the  weighting  function 
space,  we  select  a  scalar  k  and  label  modes  with  wavenumbers  k  such  that  |/cj  <  k  as  the  low 
wavenumber  or  the  coarse  scale  modes,  and  the  remaining  modes  as  the  high  wavenumber  or 
the  fine  scale  modes.  Then  depending  on  whether  an  ODE  in  the  numerical  approximation 
is  obtained  from  a  coarse  or  a  fine  scale  weighting  function,  we  classify  it  as  a  coarse  or  a 
fine  scale  equation. 

In  several  popular  methods  (such  as  the  vanishing  viscosity  method  [38])  that  guarantee 
the  convergence  of  the  numerical  solution  to  the  unique  entropy  solution,  the  numerical 
viscosity  is  applied  to  both  the  coarse  and  the  fine  scale  equations.  On  the  other  hand,  in 
the  vanishing  spectral  viscosity  method  proposed  by  Tadmor  [39],  the  viscosity  is  applied 
only  to  the  fine  scale  equations.  As  a  result,  this  method  retains  the  spectral  accuracy  of  the 
coarse  or  the  large  scale  modes  while  guaranteeing  convergence  to  the  entropy  solution.  It  is 
interesting  to  note  that  in  the  context  of  the  large  eddy  simulation  (LES)  of  incompressible 
turbulent  flows,  the  multiscale  method  of  Hughes  et  al.  [12, 14],  also  involves  applying  a 
numerical  viscosity  only  to  the  fine  scale  equations. 

Motivated  by  the  class  of  methods  where  the  viscosity  appears  only  in  the  fine  scale 
equations,  we  propose  a  method  where  different  numerical  viscosities  appear  in  the  large  and 
the  small  scale  equations.  In  addition,  in  contrast  to  the  methods  described  above,  these 
viscosities  are  not  determined  a-priori,  instead  they  are  calculated  as  part  of  the  solution 
(dynamically).  The  equations  that  are  used  to  determine  the  viscosities  are  derived  from  the 
condition  that  the  resulting  numerical  method  be  optimal  in  a  certain  user-defined  sense. 
We  dub  this  method  the  dynamic  multiscale  viscosity  method. 

We  remark  that  the  equation  used  to  dynamically  determine  the  viscosities,  is  in  effect 
the  variational  counterpart  of  the  Germano  identity.  This  identity  has  found  widespread 
use  in  determining  model  parameters  in  the  LES  of  turbulent  flows  [15].  Recently,  we 
have  demonstrated  how  it  may  be  used  as  a  tool  for  determining  unknown  parameters  in  a 
numerical  method  aimed  at  solving  an  abstract  partial  differential  equation  [16].  The  work 
presented  in  this  chapter  is  an  application  of  this  methodology  to  the  spectral  approximation 
of  non-linear  conservation  laws.  In  particular  we  use  it  to  develop  the  dynamic  multiscale 
method  for  a  generic  non-linear  conservation  law  and  then  apply  it  to  the  model  case  of  one¬ 
dimensional  Burgers  equation  to  study  its  properties.  We  find  that  the  dynamic  multiscale 
method  outperforms  the  vanishing  spectral  viscosity  method. 

An  outline  of  the  remainder  of  this  chapter  is  as  follows:  In  Section  2,  we  present  the 
equations  for  a  generic  non-linear  conservation  law.  In  Section  3,  we  introduce  its  Fourier- 
Galerkin  approximation.  In  Section  4,  we  introduce  the  multiscale  viscosity  method.  In 
Section  5,  we  derive  the  necessary  conditions  that  ensure  the  multiscale  viscosity  method  is 
optimal  in  a  user-defined  sense.  We  employ  these  conditions  in  Section  6  to  derive  expressions 
for  the  multiscale  viscosities.  This  completes  the  description  of  the  dynamic  multiscale 
method.  In  Section  7,  we  apply  the  proposed  method  to  the  one-dimensional  Burgers 
equation,  and  compare  the  results  with  the  Fourier-Galerkin  and  the  spectral  vanishing 
viscosity  methods.  We  end  with  concluding  remarks  in  Section  8. 


3.2  Problem  Statement 


We  represent  a  generic  non-linear  conservation  law  with  the  following  quasi-linear  partial 
differential  equation  in  the  space-time  domain  Q  =  flx]0,  T[,  where  Cl  =  [0, 2ir\  is  the 
spatial  domain,  and  ]0,  T[  denotes  the  time  period  of  interest:  Given  /,  G,  D  and  Uq,  find 
u  :  Q  — i ►  Rn,  such  that 

uit  +  G(u)>x  -  Duxx  +  f  =  0,  in  Q  (3.1) 

u(x,  0)  =  Uq(x),  in  Cl.  (3.2) 

In  the  equations  above,  u<t  is  the  time  derivative  of  u,  G  :  Rn  — ►  IR”  is  a  non-linear  vector 
function  of  u,  D  :  Kn  — ►  Rn  is  an  n  x  n  positive  semi-definite  matrix  of  viscosities  given 
by  D  =  diag{t'i,  •  •  •  ,  un},  and  uQ  is  the  initial  condition.  The  form  of  equations  (3. 1-3.2)  is 
representative  of  a  large  class  of  physical  phenomena  that  includes  the  dynamics  of  gasses, 
models  of  traffic  flow,  non-linear  water  waves,  and  processes  described  by  Burgers  equation. 
We  consider  periodic  boundary  conditions  for  u  expressed  as 

u(2n,  t)  =  u(0,t),  t  €]0,T[.  (3.3) 

For  any  v  :  Q  — >  Rn,  we  introduce  a  Fourier-series  representation  defined  as 

p(«)v=  v(k,t)eikx,  (3.4) 

0<|fc|<a 

where  k  is  the  wavenumber,  and  u  are  the  Fourier  coefficients  given  by 

v(k,t)  = f  v(x,  t)e~lkx  dx.  (3.5) 

Jn 

It  is  easily  verified  that  the  operator  P1^  commutes  with  spatial  and  temporal  differentiation 
and  that 


p(«)p(/3)  _  p(/3)p(“)  _  p(min(a,0))^  (3-6) 

We  will  make  use  of  this  property  in  deriving  the  consistency  conditions  for  the  optimal 
numerical  method  in  Section  5. 

We  are  interested  in  spectral  approximation  of  (3.1-3. 2)  for  small  viscosities  (|£)|  1) 

and  in  the  limit  of  vanishing  viscosities  (|23>|  — 0).  For  small  viscosities  the  solution  to  these 
equations  is  known  to  exhibit  sharp  variations  in  space  and  time  called  shocks.  For  D  =  0 
the  solution  becomes  discontinuous  and  multivalued.  To  ensure  uniqueness,  (3. 1-3.2)  must 
be  supplemented  with  conditions  that  relate  jumps  in  conserved  quantities  across  a  shock 
and  an  entropy  production  inequality.  Another  mechanism  to  arrive  at  the  same  physically 
relevant  solution  in  this  limit  is  to  construct  a  solution  with  finite  viscosity  and  then  consider 
the  limit  |G|  — >  0. 
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3.3  Fourier- Galer kin  Approximation 

The  Fourier-Galerkin  approximation  of  (3. 1-3.2)  is  a  function 

u(N)  =  u(Kt)eikx  (3.7) 

0<|/c|<N 

such  that 

P(W)  +  G(uN),x  -  Dunxx  +  /]  =0,  inQ  (3.8) 

pWuW(x,0)  =  P{N)u0(x),  in  a  (3.9) 

Noting  that  P^  commutes  with  spatial  and  temporal  differentiation,  and  that 
pMw(W)  _  U(N)  ^  (3  8-3. 9)  may  be  simplified  as 

uf]  +  (pw[G(uw)]j  -Dug  +  P(w)/  =  0,  in  Q  (3.10) 

u(N)(x,0)  =  P(N)u0(x),  inQ.  (3.11) 

Using  (3.7)  in  (3.10-3.11)  and  invoking  the  orthogonality  of  elkx  in  Q,  we  have 

uNtt  +  ikG{uw)  +  k2Du{N)  +  f  =  0,  0  <  \k\  <  N,  in  ]0,T[  (3.12) 

uw(k,  0)  =  M0(fc),  0<\k\<N,  (3.13) 


Note  that  in  (3.12),  for  notational  convenience,  we  have  omitted  the  explicit  dependence  of 
the  Fourier  coefficients  (the  hat  terms)  on  k  and  t. 

Equations  (3.10-3.11)  and  (3.12-3.13)  are  both  expressions  of  the  Fourier-Galerkin 
approximation  to  the  original  partial  differential  equation.  We  wish  to  solve  these  equations 
in  the  limit  of  small  or  vanishing  viscosities  when  the  mesh  size  denoted  by  -n/N  is  much 
larger  than  the  shock  width.  In  this  case,  the  Fourier-Galerkin  is  known  to  produce  large 
spurious  oscillations  and  become  unstable.  Further,  in  the  limit  D  =  0,  even  with  mesh 
refinement,  that  is  N  — >  oo,  the  Fourier-Galerkin  solution  is  known  not  to  converge  to 
“physical”  entropy  solution.  In  order  to  address  these  issues,  in  the  following  section  we 
propose  a  multiscale  method  based  on  adding  numerical  viscosities  to  the  Fourier-Galerkin 
approximation. 

3.4  Multiscale  Viscosity  Method 

We  augment  the  Fourier-Galerkin  method  with  multiscale  viscosities.  The  choice  of  using 
two  distinct  viscosities,  one  for  the  coarse  scale  equations,  and  another  for  the  fine  scale 
equations  is  motivated  by  the  earlier  work  of  several  researchers  [12, 14,39,40].  The  resulting 
method  is  given  by 

—  / 

u(f)  +  (pW[G(uW)])  ~(D  +  +  ^(1  -  P(oN)))[u^]  +  P{N)f  =  0, 

in  Q,  (3.14) 
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where  a  €]0, 1[  is  a  real  number,  I  is  the  identity  operator,  and  D  =  diag{P1?  •  •  •  ,  Dn}  and 
D  =  diag{zd,  •  •  •  ,  Pn}  are  matrices  of  numerical  viscosities.  The  initial  condition  remains 
unaltered  and  is  given  by  (3.11).  Note  that  two  distinct  numerical  viscosities,  given  by  D/N 
and  D/N ,  appear  in  the  equations  for  the  coarse  and  fine  scales  respectively.  This  is  clearly 
seen  once  equations  for  the  Fourier  coefficients  of  are  evaluated.  That  is,  using  (3.7)  in 
(3.14), 

uN,t  +  ikG(uW)  +  k\D  +  ~)u{N)  +  f  =  0,  0  <  |fc|  <  aN,  in  ]0, T[  (3.15) 

uN,t  +  ikG(uW)  -f  k\D  +  ^)u{N)  +  f  =  0,  aN<\k\<  N,  in  ]0,  T[.  (3.16) 

Since  (3.15)  is  obtained  from  coarse  scale  weighting  functions  ( \k\  <  aN),  and  (3.16)  is 
obtained  from  fine  scale  weighting  functions  we  term  these  equations  the  coarse  and  fine 
scale  equations  respectively.  As  is  apparent  from  these  equations,  the  numerical  viscosity 
that  appears  in  the  coarse  and  fine  scale  equations  is  given  by  D/N  and  D/N  respectively. 

While  the  proposed  method  is  motivated  by  earlier  works  (see  [12,14,39,40]),  there  are 
two  crucial  differences 

1.  In  the  aforementioned  methods  the  viscosity  is  applied  only  to  the  fine  scale  equations. 
The  viscosity  in  the  coarse  scale  equations  is  zero.  In  our  method  different  but  non-zero 
viscosities  are  applied  to  both  the  coarse  and  the  fine  scale  equations. 

2.  In  the  aforementioned  methods  the  viscosity  is  determined  a-priori.  In  our  method 
the  viscosity  is  determined  as  a  part  of  the  calculation  using  an  optimality  argument. 
This  development  is  described  in  the  following  section. 

з. 5  Consistency  Conditions 

We  now  derive  a  set  of  consistency  conditions  that  are  utilized  in  the  next  section  to  compute 
an  explicit  expression  for  evaluating  the  numerical  viscosities.  These  conditions  are  motivated 
by  the  Germano  identity,  which  is  commonly  used  to  evaluate  parameters  in  subgrid  models 
for  the  large  eddy  simulation  of  turbulent  flows  [15]. 

The  main  idea  which  is  expressed  in  Theorem  5.1  below,  is  the  following.  We  assume  that 
it  is  possible  to  choose  the  viscosities  in  (3.14)  such  that  the  resulting  solution  is  optimal 
in  the  sense  that  its  Fourier  coefficients  exactly  match  the  corresponding  coefficients  of  the 
continuous  solution.  Note  that  other,  user-defined  definitions  of  an  optimal  solution  are  also 
possible.  In  addition,  we  assume  that  the  same  viscosities  also  yield  optimal  results  for  a 
coarser  discretization.  That  is  the  solution  of  (3.14)  with  N  replaced  by  M  everywhere, 
where  M  <  N,  is  also  optimal  in  the  manner  described  above.  These  assumptions  lead  to  a 
set  of  conditions  that  must  be  satisfied  by  the  numerical  viscosities  in  order  to  yield  optimal 
results.  A  key  feature  of  these  conditions  is  that  they  do  not  involve  the  continuous  solution 

и,  and  are  expressed  entirely  in  terms  of  the  numerical  solution  uN.  Hence  they  may  be 
utilized  with  relative  ease  to  determine  the  numerical  viscosities. 


Theorem  5.1:  Let  and  be  solutions  of  the  multiscale  viscosity  method  with 
modes  up  to  N  and  M  respectively,  with  M  <  N.  If 


U(N)  =  p  {N)u 
U{M)  =  p  (M)u 

where  u  is  the  solution  of  the  continuous  problem,  then 


(3.17) 

(3.18) 


p(orW) 

N 


)  +  D(- 


_ jp(aM) 

~M 


))[(p(M)uW)ix]  = 
^p(M)[G(p(M)u(W))  _  G(U(W))]) 


,X 


in  ($3.19 ) 


Proof:  vSMS>  satisfies  (3.14)  with  N  replaced  by  M  everywhere.  Further,  from  (3.18) 

U(M)  =  p(M)u 

=  FiM)Fwu  (from  (3.6),  and  since  M  <  N) 

=  FiM)u{N)  (from  (3.17)).  (3.20) 

Using  (3.20)  in  (3.14)  written  with  N  replaced  by  M,  and  rearranging  terms  so  as  to  retain 
only  the  model  term  on  the  left-hand  side,  we  arrive  at 

(§P-  +  ^(1  -  P(“M)) )  =  P(MWtw)  +  (p(m)[G(P(m)w(w))]) 

-£>(P(m)uw),xx  +  P(M)/,  in  Q.(3.21) 


By  applying  P^M*  to  (3.14),  assuming  that  P^M^  and  spatial  differentiation  commute, 
using  property  (3.17),  and  rearranging  terms  so  as  to  retain  only  the  model  term  on  the  left 
hand  side,  we  conclude  that 

(J^F(aN)  +  ^(1  -  F{aN))^j  [(P(M)uw),xx]  =  P(M)u(fW)  +  (p (M)[G(u(JV))]) 

-DF{m)u^  +  P in  Q.  (3.22) 

Subtracting  (3.22)  from  (3.21)  we  have  the  desired  result  (viz.  (3.19))  ■ 

Remark:  The  Fourier  representation  of  (3.19)  is  given  by 


-K  < 


'  (§  -  .0  <  |fc|  <  aM 

D  D 


(rL_^,)(&<">)  ,aM<\k\<(3  ^=iA:(G(p(M)WW)-G(^)))in]0,r[(3.23) 


where 


j3  =  min(aiV,  M ). 


(3.24) 
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3.6  Evaluation  of  the  Viscosity  Parameters 

Equation  (3.19)  (or  (3.23))  represents  as  many  relations  as  there  are  modes  for  which 
|  A:  |  <  M.  We  wish  to  evaluate  only  the  2 n  parameters  that  appear  in  D  and  D  using 
these  relations.  One  mechanism  of  reducing  (3.19)  to  2 n  equations  is  to  equate  the  L2  inner 
product  of  its  residual  with  the  linearly  independent  functions 

p(aM)uWj  j  =  (3.25) 


and 


(I  _  P(«M))U(JV))  j  =  lf . . .  >  n  (3  26) 

to  zero.  This  procedure  leads  to  2n  scalar  equations  that  are  most  conveniently  expressed 
in  terms  of  the  Fourier  coefficients  of  For  j  =  1,  •  •  •  ,n  they  are  given  by 

53  (3.27) 

0<\k\<ccM 


e  «r«(  a,-[P(M)uW]  -  Gjtw(W)])(3.28) 

aM<\k\<M 

where  the  repeated  j  indices  do  not  imply  a  summation.  These  expressions  for  evaluating 
the  viscosity  parameters  are  functions  of  the  solution  itself.  Thus  the  closed  system  of 
equations  (3.14),  (3.11)  and  (3.27-3.28)  comprises  a  numerical  method  with  non-linear  (in 
Wat),  multiscale  viscosities.  In  implementing  this  method  the  viscosities  are  evaluated  based 
on  the  solution  obtained  from  the  previous  time-step.  However  once  the  numerical  viscosity 
is  determined  this  term  is  treated  implicitly. 


-<!-£)  E  W  - 


0<|fc|<aM 


4-|>  E  *2i“f’i2- 


aM<\k\<0 


<§-$  ,e  *2i“ri2  = 


0<\ k\<M 


3.7  Numerical  Example 

As  an  example  we  apply  the  proposed  method  to  Burgers  equation  in  one  dimension.  In  this 
case,  in  (3. 1-3. 2),  n  =  1,  /  =  0,  G(u)  —  u2/ 2,  D  —  iq  =  5  x  10-5,  and  u0  =  —  sinx. 

The  solution  to  this  problem  evolves  in  two  distinct  phases.  In  the  first  phase  ( t  <  7t/2)  , 
the  smooth  sine  curve  steepens.  In  wavenumber  space,  this  corresponds  to  transfer  of  energy 
from  the  k  —  1  mode  to  higher  wavenumbers.  This  phase  culminates  in  the  formation  of  a 
shock  (or  an  inverted  N-wave)  at  t  «  7r,  whose  approximate  width  is  l  =  1.6  x  10-4  units. 
All  the  dissipation  in  the  system  is  concentrated  near  the  shock.  In  the  wavenumber  space, 
the  formation  of  the  shock  corresponds  to  a  1  fk  spectrum  that  extends  to  the  dissipation 
wavenumber,  where  it  steepens.  In  the  second  phase,  the  magnitude  of  the  N-wave  reduces 
as  1/(1  +  t)  as  the  strength  of  the  shock  weakens.  In  wavenumber  space,  this  corresponds 
to  the  lowering  of  the  entire  spectra  at  the  rate  of  1/f.  These  stages  of  the  solution  are 
presented  in  Figures  3.1  and  3.2,  where  we  have  shown  a  well-resolved  numerical  solution  of 
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the  problem.  Note  that  E(k)  =  represents  the  magnitude  of  the  Fourier  coefficients  of 
the  solution. 

In  order  to  asses  the  performance  of  the  proposed  method  we  consider  the  following 
numerical  solutions 

1.  A  Fourier-Galerkin  solution  obtained  by  solving  (3.10-3.11)  or  (3.12-3.13)  on  a  fine  grid 
in  which  the  shock  is  resolved.  Following  terminology  used  in  turbulence  modeling  we 
refer  to  this  benchmark  solution  as  the  direct  numerical  solution  (DNS).  The  number 
of  modes  used  for  computing  the  DNS  solution  is  N  =  65, 536.  This  corresponds  to  a 
mesh  size  h  =  4.79  x  10“5,  which  is  smaller  than  the  shock  width  l  =  1.6  x  10~4.  In 
Figures  3.1  and  3.2  we  have  plotted  this  solution  in  physical  and  wavenumber  spaces 
at  various  instances  during  the  interval  ]0,5[. 

2.  A  Fourier-Galerkin  solution  obtained  by  solving  (3.10-3.11)  or  (3.12-3.13)  on  a  coarse 
grid  with  N  =  64.  In  this  case  the  finest  resolved  scale  (h  =  4.9  x  10~2)  is  much  coarser 
than  the  scale  at  which  dissipation  occurs  ( l  —  1.6  x  10“4). 

3.  A  vanishing  spectral  viscosity  solution  on  a  coarse  grid  with  N  =  64.  This  method  is 
represented  by  (3.14),  where  D  =  Pi  =  0,  and  D  is  non-zero.  In  particular  we  choose 
ct  =  1/2  and  D  —  id  =  0.25a:.  This  choice  for  v\  is  based  on  the  guideline  provided 
in  [39]. 

4.  A  dynamic  multiscale  viscosity  solution  on  a  coarse  grid  with  N  =  64.  This  method 
is  given  by  (3.14),  where  a  =  1/2  and  the  viscosity  coefficients  are  determined  using 
(3.27-3.28).  In  calculating  these  coefficients  the  solution  from  the  previous  time-step 
is  used. 

3.7.1  Comparison 

The  viscosity  parameters  of  the  dynamic  multiscale  method  are  chosen  such  that  the  method 
satisfies  conditions  that  are  necessary  to  ensure  that  the  resulting  numerical  solution  has  the 
same  Fourier  coefficients  as  the  continuous  solution.  This  stipulation  may  be  interpreted 
as  a  criterion  used  to  design  the  numerical  method.  In  the  following  comparison  we  assess 
how  close  the  method  comes  to  achieving  this  criterion.  We  also  assess  its  performance  in 
relation  to  other  methods. 

In  Figures  3.3-3.6,  we  have  plotted  E(k)  for  the  three  numerical  methods  and  the 
truncated  DNS  solution  at  four  distinct  times.  The  DNS  serves  as  a  benchmark  solution.  At 
t  =  0.5  we  observe  that  there  is  little  difference  in  the  numerical  methods  and  the  DNS.  At 
t  =  1.5,  when  the  shock  is  about  to  form  differences  appear.  A  significant  pile-up  of  energy 
near  the  cut-off  wavenumber  is  observed  in  the  coarse  Fourier-Galerkin  solution.  For  the 
vanishing  spectral  viscosity  solution,  this  pile-up  is  reduced.  However  the  solution  is  seen  to 
oscillate  about  the  DNS  at  wavenumbers  close  to  the  separation  between  the  coarse  and  fine 
scales  ( k  =  32).  The  multiscale  solution  has  much  smaller  oscillations  at  these  wavenumbers 
and  only  slightly  underestimates  the  spectrum  at  higher  wavenumbers.  At  time  t  =  2,  we 
observe  that  the  pile-up  at  high  wavenumbers  in  the  Fourier-Galerkin  solution  has  polluted 
the  results  at  lower  wavenumbers.  The  vanishing  spectral  viscosity  solution  continues  to  be 
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accurate  at  the  lower  wavenumbers,  however  the  oscillations  close  to  k  =  32  appear  to  have 
increased.  The  dynamic  multiscale  viscosity  solution  is  accurate  at  the  lower  wavenumbers. 
The  oscillations  near  k  =  32  persist,  however  they  are  less  pronounced  than  those  for  the 
vanishing  spectral  viscosity  solution.  At  £  —  5  we  observe  that  Galerkin  solution  is  completely 
inaccurate,  the  oscillations  in  the  vanishing  spectral  viscosity  solution  have  increased  and 
propagated  to  lower  wavenumbers,  where  as  the  dynamic  multiscale  solution  has  retained 
its  accuracy. 

In  Figure  3.7,  we  have  plotted  the  viscosity  parameters  Pi/N  and  V\jN  for  the  vanishing 
spectral  viscosity  method  and  the  dynamic  multiscale  method  as  a  function  of  time.  The 
coarse  scale  parameter  for  the  vanishing  spectral  viscosity  method  is  zero  and  is  not  shown, 
whereas  the  fine  scale  parameter  which  is  non-zero  and  constant  is  shown.  For  the  dynamic 
multiscale  method  we  observe  that  both  the  coarse  and  fine  scale  parameters  are  zero  till 
£  ~  1.  This  represents  the  time  it  takes  for  the  energy  to  spill  out  to  the  wavenumbers  near 
the  cutoff  wavenumber.  Thus  the  numerical  method  (correctly)  imposes  no  viscosity  till  this 
time.  Thereafter,  the  fine  scale  parameter  is  seen  to  rise,  and  after  a  while  the  coarse  scale 
parameter  follows  suite.  A  couple  of  observations  are  noteworthy:  1)  Unlike  the  vanishing 
spectral  viscosity  method,  the  viscosity  in  the  coarse  scales  in  the  dynamic  multiscale  method 
is  not  zero,  thus  the  method  is  qualitatively  different;  2)  The  viscosity  in  the  coarse  and  the 
fine  scales  is  active  for  different  periods  of  time,  and  also  has  different  values.  In  particular, 
the  fine  scale  viscosity  is  about  5/3  of  the  coarse  scale  value. 

Next  we  compare  the  accuracy  with  which  the  numerical  methods  predict  the  decay  of 
kinetic  energy  with  time.  We  define  the  relative  error  in  the  resolved  kinetic  energy  as  follows 


ek e(£)  — 


W64)N2 


(3.29) 


In  Figure  3.8  we  have  plotted  ei«,(£)  for  all  the  numerical  methods.  We  observe  that  at 
£  ~  1,  that  is  when  the  spectra  begins  to  spill  beyond  the  numerical  cut-off  wavenumber, 
the  Fourier-Galerkin  solution  overestimates  the  kinetic  energy.  By  £  w  2.5  its  kinetic 
energy  is  about  two  times  the  actual  value.  Both  the  vanishing  spectral  viscosity  and  the 
dynamic  multiscale  methods  are  much  more  accurate.  In  Figure  3.9  we  have  compared  the 
performance  of  these  two  methods.  We  observe  that  the  vanishing  spectral  viscosity  solution 
overestimates  the  kinetic  energy  and  that  the  error  is  seen  to  increase  with  time.  At  £  =  5  the 
total  error  is  about  2.7%.  The  multiscale  solution  underestimates  the  kinetic  energy  however 
the  error  is  less  (about  0.2%  at  £  =  5)  and  remarkably  it  does  not  appear  to  increase  with 
time. 

One  of  the  attractive  features  of  the  vanishing  spectral  viscosity  method  is  that  in  the  limit 
v  — »  0,  it  converges  to  the  unique  solution  that  satisfies  the  entropy  production  inequality, 
while  retaining  spectral  accuracy  in  the  coarse  modes.  This  leads  us  to  consider  that  while 
the  dynamic  multiscale  method  may  be  more  accurate  in  predicting  quantities  such  as  the 
overall  spectrum  and  the  resolved  kinetic  energy,  the  vanishing  spectral  viscosity  method 
may  be  more  accurate  in  capturing  the  evolution  of  the  coarse  modes  since  it  imposes  no 
additional  viscosity  on  these  modes.  In  order  to  verify  this  in  Figure  3.10,  we  have  plotted 
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(3.30) 


the  error  in  the  k  =  1  mode,  scaled  by  the  exact  value,  as  a  function  of  time.  That  is 

u^64)(l,t)  —  u(l,t) 

We  consider  the  DNS  solution  to  be  the  exact  value.  We  observe  that  at  t  «  7r  the  error  in 
the  coarse  Fourier-Galerkin  solution  rises  steeply,  whereas  the  error  in  the  other  two  methods 
is  smaller.  In  Figure  3.11,  we  exclude  the  Fourier-Galerkin  solution.  We  observe  that  error 
the  vanishing  spectral  viscosity  solution  rises  steadily  beyond  t«l,  At  t  =  5  the  error  is 
about  1%.  The  error  in  the  dynamic  multiscale  method  is  much  smaller  and  appears  not 
to  increase  with  time.  At  t  =  5  it  is  about  0.04%.  We  have  observed  similar  behavior  for 
other  coarse  modes  (k  =  1,  •  •  •  ,  10).  Thus  contrary  to  what  might  be  expected,  we  observe 
that  dynamic  method  (with  a  non- zero  viscosity  in  the  coarse  modes)  is  more  accurate  than 
the  vanishing  spectral  viscosity  method  in  predicting  the  evolution  of  the  coarse  modes.  We 
attribute  this  observation  to  the  hypothesis  that  the  ideal  model,  which  would  replicate  the 
effects  of  the  missing  scales  on  the  retained  scales  exactly,  may  possess  a  non-zero  viscosity  at 
small  wavenumbers.  We  are  presently  verifying  this  hypothesis  analytically  and  numerically. 


<h(t) 


3.8  Conclusions 

We  have  proposed  a  new  dynamic  multiscale  viscosity  method  for  the  spectral  approximation 
of  conservation  laws  in  the  limit  of  small  or  vanishing  viscosities.  Within  this  method 
the  numerical  approximation  is  split  into  coarse  and  fine  scales,  and  likewise  the  projected 
spectral  equations  are  also  split  into  coarse  and  fine  scale  equations.  Thereafter  different 
numerical  viscosities  are  applied  in  the  coarse  and  fine  scale  equations.  These  viscosities  are 
determined  using  a  condition  that  must  be  satisfied  if  the  resulting  numerical  solution  is  to 
be  optimal  in  a  user-defined  sense. 

We  have  applied  this  method  to  the  one-dimensional  Burgers  equation.  We  have 
compared  the  resulting  solution  with  the  Fourier-Galerkin  solution  and  the  vanishing  spectral 
viscosity  solution  computed  using  the  same  number  of  modes.  As  a  benchmark  we  have  used 
a  well-resolved  Fourier-Galerkin  solution.  In  all  comparisons  we  have  found  that  the  dynamic 
multiscale  solution  is  the  most  accurate.  In  addition  we  have  observed  that  the  relative  errors 
in  this  solution  appear  not  to  grow  in  time. 
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Figure  3.1:  Well  resolved  numerical  solution  (D 

t  =  0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5  and  5.0  (arrow  indicates  increasing  time) 
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Figure  3.2: 
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Figure  3.5:  E(k)  for  the  solution  of  the  numerical  methods  and  the  DNS  at  t  —  2.0. 
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Figure  3.7:  Numerical  viscosities  for  the  spectral  vanishing  viscosity  method  and  the  dynamic 
multiscale  method  as  a  function  of  time. 


Figure  3.8:  Scaled  error  in  the  resolved  kinetic  energy  of  the  numerical  methods  as  a  function 
of  time. 
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Chapter  4 


Variational  Germano  Identity  Applied 
to  Large  Eddy  Simulation 

4. 1  Introduction 

The  typical  formalism  for  performing  large  eddy  simulation  (LES)  of  turbulent  flows  involves 
spatially  filtering  the  strong  form  of  the  Navier-Stokes  equations.  The  width  of  the  filter  is 
chosen  in  accordance  with  the  resolution  of  the  computational  grid  on  which  the  numerical 
solution  is  sought.  The  filtering  operation  yields  a  system  of  equations  for  the  filtered 
variables  in  which  certain  terms  that  involve  the  unfiltered  variables  appear.  These  terms 
are  replaced  by  models,  which  are  functionals  of  the  filtered  variables  only,  and  a  closed 
system  for  the  filtered  variables  is  obtained. 

In  [14],  a  different  formalism  for  performing  LES  was  introduced.  In  this  new  approach, 
instead  of  spatially  filtering  the  strong  form  of  the  Navier-Stokes  equations,  a  weak  or  a 
variational  form  was  used  as  a  starting  point.  It  was  recognized  that  the  approximation  of 
the  weak  form  by  the  Galerkin  method  involved  approximating  infinite-dimensional  function 
spaces  by  their  finite-dimensional  counterparts.  It  was  argued  that  for  typical  turbulent  flows, 
which  tend  to  have  a  large  range  of  energetic  structures  of  different  sizes,  finite  dimensional 
spaces  could  not  faithfully  represent  all  the  features  of  the  exact  solution.  In  the  few  select 
cases  where  this  would  be  feasible,  the  Galerkin  solution  would  correspond  to  the  direct 
numerical  simulation  of  the  problem.  On  the  other  hand,  in  most  cases  of  practical  interest, 
LES  models  that  represent  the  effect  of  the  missing  scales  will  be  required.  This  formalism 
was  then  used  to  develop  a  multiscale  model  (in  reality,  a  two-scale  model)  for  the  resolved 
scales.  This  model  consisted  of  a  Smagorinsky-type  term  constructed  from  the  resolved, 
fine-scale  velocity  components  and  applied  to  the  fine-scale  equations  and  no  model  in  the 
coarse  scale  equations.  Several  applications  of  this  model  with  encouraging  results  have  now 
been  reported  (see  for  example  [12,13,41]). 

While  the  variational  formulation  of  LES  has  lead  to  encouraging  results,  there  is  one 
aspect  of  this  formulation  that  is  still  missing.  This  is  the  development  of  dynamic  variational 
LES  models.  In  the  context  of  filter-based  LES  modeling,  dynamic  LES  models  rely  on  the 
Germano  identity  to  evaluate  model  parameters  as  part  of  the  simulation  (dynamically)  [15]. 
Dynamic  LES  models  have  been  particularly  successful  in  modeling  flows  where  the  statistics 
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of  turbulence  vary  in  space  and  time,  such  as  flows  with  intermittent,  transitional  or  decaying 
turbulence  and  flows  with  boundary  layers.  In  this  manuscript,  we  derive  the  variational 
counterpart  of  the  Germano  identity  and  demonstrate  how  it  may  be  used  to  determine  LES 
model  parameters  in  a  variational  context. 

Our  derivation  of  the  variational  counterpart  of  the  Germano  identity  is  similar  to  its 
derivation  in  the  filtered  case.  There  are  however  important  differences  in  the  final  result. 
These  differences  emanate  from  differences  in  the  definition  of  an  ideal  LES  model  in  the 
variational  and  the  filter-based  formulations  and  are  highlighted  in  our  development.  In  fact, 
we  have  discovered  that  the  variational  formulation  of  the  Germano  identity  has  applications 
beyond  the  realm  of  large  eddy  simulation.  In  [16]  we  have  demonstrated  how  it  may  be  used 
to  determine  parameters  of  a  generic  numerical  method.  So  far  in  addition  to  its  application 
to  turbulence  (described  in  this  chapter),  we  have  applied  it  to  the  determine  parameters  in 
a  residual-based  finite  element  method  solution  of  the  linear  advection-diffusion  problem  [16] 
and  to  determine  the  multiscale  viscosity  parameters  in  the  spectral  approximation  of  Burgers 
equation  [17]. 

The  layout  of  the  remainder  of  this  chapter  is  as  follows:  In  the  following  section  we 
introduce  the  incompressible  Navier  Stokes  equations.  In  Section  3,  we  review  the  filter- 
based  Germano  identity.  Thereafter  in  Section  4,  we  introduced  the  variational  formulation 
of  LES  and  also  derive  the  variational  counterpart  of  the  Germano  identity.  In  Section  5,  we 
apply  it  to  the  specific  problem  of  determining  the  Smagorinsky  parameter  in  the  spectral 
approximation  of  the  decay  of  homogeneous  isotropic  turbulence.  In  Section  6,  we  present 
numerical  results  and  end  with  conclusions  in  Section  7. 


4.2  Incompressible  Navier  Stokes  Equations 

We  consider  the  motion  of  an  incompressible  fluid  in  Q  =  flx]Ti,  T2[,  where  Cl  is  the  spatial 
domain  and  ]Ti,T2[  is  the  time  period  of  interest.  The  strong  form  of  the  problem  is  given 
by 

uyt  +  V  •  (it  <8>  u)  +  Vp  —  =  /,  in  Q  (4.1) 

V  •  u  =  0,  in  Q.  (4.2) 

In  the  above  equations  U  =  [u,p]T,  where  u  is  the  fluid  velocity  field  and  p  is  the  pressure. 

In  addition,  v  is  the  kinematic  viscosity  and  the  symbol  <g>  denotes  the  outer  product  of  two 

vectors.  The  boundary  of  Cl  is  denoted  by  dCl.  For  simplicity  the  boundary  condition  on  U 
is  assumed  to  be  either  periodic  or  homogeneous.  The  initial,  divergence-free  velocity  field 
is  given  by 

u(x,  0)  =  uq(x).  (4.3) 


4.3  Filtered  Form  of  the  Germano  identity 

In  this  section  we  provide  a  summary  of  the  derivation  the  filtered  form  of  the  Germano 
identity  (for  details  see  [15]).  This  will  be  contrasted  with  the  its  variational  counterpart  in 
the  following  section. 
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We  apply  a  spatial  filter  to  (4.1)  and  (4.2)  to  arrive  at 

(F hu)<t  +  V  •  F *(u  <g>  u)  +  V(F hp)  -  uV2(Fhu)  =  F hf,  in  Q  (4.4) 

V  •  Fhu  =  0,  inQ  (4.5) 

where  F,lw  is  the  filtered  version  of  u  and  ¥h  is  a  filter  of  width  h  and  we  have  assumed 
that  filtering  and  spatial  differentiation  commute.  Note  that  (4.4)  and  (4.5)  are  not  closed 
in  Fhw  as  (4.4)  still  involves  the  unfiltered  velocity  u  in  the  quadratic  term.  Thus  (4.4)  may 
be  replaced  by 

(Fhu),t  +  V  •  (Fftw  0  Fftw)  +  V(FV)  -  isV2(Fhu)  -  V  •  m(¥hu ;  h,  c)  =  ¥hf,  in 

where  m(Fhu ;  h,  c )  is  a  model  term  that  depends  only  on  the  filtered  velocity  field  ¥hu  (and 

not  u ),  the  filter  width  h  and  a  vector  of  parameters  c.  For  a  fixed  choice  of  c,  (4.6)  and 

(4.5)  represent  a  closed  system  in  F hu  and  ¥hp. 

From  (4.6)  and  (4.4)  we  conclude  that  the  ideal  model  term  is  given  by 

m(¥hu]h,c)  =  (¥ku)  <S>  (¥hu)  -  ¥h(u  0  u)  in  Q,  (4.7) 

as  this  choice  allows  the  filtered  velocity  and  pressure  fields  [F hu,  F/lp]  to  be  the  solution  of 

(4.6)  and  (4.5). 

To  derive  the  filtered  form  of  the  Germano  identity  we  consider  a  coarser  spatial  filter 
FH ,  with  H  >  h.  Using  the  argument  made  in  the  previous  paragraph,  we  conclude  that  at 
this  scale  the  ideal  model  is  given  by 

m(Fffu;  H,  c )  =  (FHu)  ®  (FHu)  —  F^-u  <g>  u )  in  Q.  (4.8) 

Note  that  we  have  assumed  that  the  functional  dependence  of  the  model  on  the  filtered 
solution  and  the  filter  width  as  well  as  the  parameters  c  is  unchanged.  By  subtracting  (4.7) 
filtered  at  the  H  scale  from  (4.8),  and  assuming  that  ¥H  =  F/fF/l,  we  arrive  at  the  Germano 
identity,  viz, 

m(FHw;  H,  c)  -  FH(m(FV  h,  c))  =  (F Hu)  0  (F*u)  -  F*((F hu)  0  {¥hu))  in  Q(4.9) 

This  expression  only  involves  the  filtered  velocity  fields  and  does  not  contain  any  unfiltered 
quantity.  It  was  first  introduced  in  [15]  and  since  then  has  found  wide  application  in 
determining  parameters  (denoted  by  c  in  our  case)  in  LES  models.  Note  that  strictly 
speaking  this  relation  is  valid  only  for  the  exact  filtered  field,  which  is  also  unknown  during 
a  numerical  simulation.  However,  in  typical  applications  the  numerical  approximation  of  the 
filtered  field  is  used  in  this  relation. 

The  derivation  of  the  Germano  identity  described  above  is  simplified  by  making  the 
assumptions  that  spatial  differentiation  and  filtering  operations  commute  and  that  F;/F/l  = 
¥H .  Neither  of  these  assumptions  are  typically  valid  for  flows  that  are  inhomogeneous 
(wall  bounded  flows  for  example).  For  such  flows,  the  assumption  of  commutativity  of 
differentiation  and  filtering  introduces  errors  in  the  modeled  equations  which  can  be  quite 
large  (see  for  example  [42]  and  [43]).  We  are  not  aware  of  a  filtered  form  of  the  Germano 
identity  that  accounts  for  these  errors.  However  filters  that  minimize  these  errors  have  been 
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proposed  (see  for  example  [44]).  T/he  assumption  ¥n¥h  =  ¥H  may  be  avoided  altogether  by 
utilizing  the  composite  filter  in  (4.8)  instead  FH.  Note  that  the  width  of  this  composite 
filter  (denoted  by  H)  may  not  equal  to  H  and  is  determined  separately.  In  this  case,  (4.9) 
is  replaced  by 

m(FH¥hu]  H,c)-Fh (m(Fhu-,  h,c))  =  {FHF hu)  ®  (FHFhu) 

-Fh({¥ hu)  <g>  (F *«))  in  Q.  (4.10) 

4.4  Variational  Germano  identity 

We  now  derive  the  variational  form  of  the  Germano  identity.  Note  that  even  though  our 
derivation  mimics  the  derivation  of  the  filtered  version  described  in  the  previous  section,  it 
is  not  restricted  to  homogeneous  flows  or  uniform  meshes.  At  every  step  we  will  point  out 
the  salient  differences  between  the  filtered  and  the  variational  approaches. 

4.4.1  Variational  Formulation  of  LES 

We  begin  with  a  weak  or  a  variational  formulation  of  the  Navier  Stokes  equations  (see  [29] 
for  example)  given  by:  Find  U(-,t)  =  [u(-,f),p(-,t)]T  G  V,  such  that 

B(W,U)  =  (W,F),  VW  =  [w,q}T  eV,Vte}TuT2[,  (4.11) 

where  the  semi-linear  form  B(-,  -)  and  is  defined  as 

B(W,  U)  =  (w,  urt)  —  (Viw,  u  <S>  u)  +  ( w ,  Vp)  -I-  2v(Vsw,  Vsu)  —  ( Vq ,  u),  (4.12) 

and  where 


(W,F)  =  (w,f).  (4.13) 

In  the  equations  above  (•,  •)  is  used  to  denote  the  L2  inner  product  in  and  Vs u  = 

i(Vu  +  Vur). 

The  infinite  dimensional  function  space  V  is  defined  as 

V  =  {V  =  [v,r]r|u  G  Hl(fl);r  G  L2(Q)  ],  (4.14) 

where  H 1  (12)  denotes  the  Sobolev  space  of  vector- valued  functions  that  are  square-integrable 
and  whose  derivatives  are  also  square-integrable.  This  implies  that  the  energy  and  the 
enstrophy  of  the  weak  solutions  remain  finite.  Also  L2(  12)  denotes  the  space  of  scalar 
functions  that  are  square-integrable.  In  addition,  we  assume  that  the  boundary  conditions 
(though  not  explicitly  stated)  are  built  into  the  definition  of  V. 

Let  Vh  G  V  be  a  finite  dimensional  subspace.  T'he  Grulerkin  method  when  applied  to 
(4.11)  yields  the  following  equation  for  the  approximate  solution  Uh.  Find  Uh(-,t )  G  Vh, 
such  that 

B(Wh,Uh)  =  (Wh,F),  VWb  G  Vhyt  g]Ti,T2[.  (4.15) 

Remarks 
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1.  Since  Vh  is  a  subspace  of  V,  (4.11)  implies  that  the  exact  solution  U  satisfies 

B(Wh,U)  =  (Wh,F),  VWh  eVhyte]TuT2[.  (4.16) 

2.  The  Galerkin  solution  essentially  corresponds  to  a  Direct  Numerical  Simulation  (DNS) 
solution  of  the  original  problem.  If  the  space  Vh  is  incapable  of  representing  the  large 
inertial  scales  of  motion  as  well  as  the  fine  dissipative  scales,  then  this  solution  is 
inaccurate.  It  may  be  improved  upon  by  adding  a  model  term  to  the  terms  that 
appear  in  the  Galerkin  approximation. 

The  addition  of  a  model  term  to  the  Galerkin  approximation  leads  to  the  following 
variational  equation  :  Find  Uh(-,t)  €  Vh,  such  that 

B(Wh,  Uh)  +  M(Wh,  Uh\  h,  c)  =  (Wh,  F),  VW'1  €  Vh,  Vt  €]Ti,T2[.  (4.17) 

In  (4.17),  M(Wh,  Uh]  h,  c)  is  the  model  term  that  is  linear  in  Wh  and  non-linear  in  Uh.  It 
also  depends  on  the  mesh  size  h  and  a  vector  of  parameters  c.  The  equation  above  is  the 
variational  counterpart  of  (4.6).  Also  note  that  due  to  the  presence  of  the  model  term,  the 
solution  of  (4.17)  and  the  Galerkin  solution  are  distinct. 

We  now  introduce  the  concept  of  an  optimal  solution  in  the  space  Vh.  For  this  we  define 
a  restriction  operator  P/l  :  V  — >  Vh,  such  that  FhU  €  Vh  is  the  optimal  representation  of  U 
in  Vh.  For  example,  FhU  may  be  the  nodal  interpolant  or  a  suitable  projection  (L2  or  Hl  in 
our  case)  of  U  in  Vh,  in  which  case  Fh  will  be  the  interpolation  or  the  projection  operator, 
respectively.  We  would  like  to  select  our  model  term  in  (4.17)  such  that  the  solution  Uh  is 
optimal.  Thus  requiring  Uh  =  P hU  in  (4.17)  we  have 

M{Wh)  FhU ;  h,  c)  =  (Wh,  F)  -  B(Wh,FhU),  VWh  €  Vh,\/t  e]Tx,T2[.  (4.18) 

Using  (4.16)  in  (4.18)  we  arrive  at 

M(Wh,FhU;h,c)  =  B{Wh,U)- B(Wh,FhU),  VWh  eVh,Vt  e]TuT2[.  (4.19) 

This  equation  defines  the  model  required  to  generate  the  optimal  solution.  By  comparing 
the  definition  of  the  ideal  model  in  the  variational  context  (4.19)  with  its  counterpart  in  the 
filtered  approach  (4.7)  we  find: 

1.  In  the  filtered  case  only  the  non-linear  terms  in  the  Navier  Stokes  equations  contribute 
to  the  definition  of  the  ideal  model.  In  the  variational  case  every  term  in  the  semi-linear 
form  contributes. 

2.  In  the  filtered  case  it  is  assumed  that  for  the  ideal  model,  (4.7)  is  satisfied  pointwise. 
Thus  the  residual  of  (4.7)  when  multiplied  with  any  weighting  function  and  integrated 
over  Q  is  equal  to  zero.  In  contrast  to  this  (4.19)  holds  only  for  weighting  functions 
contained  in  Vh.  Thus  the  constraint  on  the  model  term  in  the  variational  case  is 
weaker.  In  particular,  in  the  variational  case  the  effect  of  the  model  on  weighting 
functions  outside  ofVh  is  irrelevant. 

3.  In  the  filtered  case  the  model  term  is  a  tensor.  In  the  variational  case  it  appears  in 
the  weak  form  just  as  a  vector  would  appear. 

These  differences  in  the  definition  of  the  ideal  model  in  the  variational  and  filtered 
approaches  are  important  and  will  lead  to  different  forms  of  the  Germano  identity. 
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4.4.2  The  variational  Germano  Identity 

Next  we  consider  another  function  space  VH  such  that  VH  C  Vh  C  V.  We  also  define  a 
restriction  operator  Pff  such  that  FH  :  V  — >  VH  and  P HU  E  VH  is  the  optimal  representation 
of  U  in  VH.  We  assume  that  the  operators  Fh  and  FH  are  related  and  hence  the  same  form 
of  the  model  term  yields  optimal  solutions  in  and  VH.  The  finite  dimensional,  modeled 
variational  equation  for  this  space  is  given  by:  Find  UH(-,t)  G  VH ,  such  that 

b(wh,uh)  +  m(wh,  uh\  h,  c)  =  ( wh,f %  vwH  e  vHyt  g]t1}t2[,  (4.20) 

where  the  same  vector  of  parameters  appears  in  (4.20)  and  (4.17).  Hence  the  model  that 
yields  the  optimal  solution  UH  =  PHU,  satisfies 

m(wh,fhu-,h,c)  =  b(wh,u)  -  b(wh,fhu),  'iwH  g  vHyt  €]T1,r2[.  (4.21) 

Returning  to  (4.19),  we  note  that  it  holds  for  all  Wh  G  V*,  and  hence  is  also  valid  for  all 
WH  G  VH .  Thus  we  have 

M(WH ,FhU\h,c)  =  B(Wh,U)  -  B(WH,FhU),  VWH  G  VH,Vt  e]TuT2[.  (4.22) 

Subtracting  (4.22)  from  (4.21)  we  arrive  at 

M(WH,FHU;H,c)-M(WH,FhU;h,c)  =  B(WH,PhU)  -  B(WH,FHU), 

VWH  G  VHyt  e]TuT2[.  (4.23) 

In  addition  if  the  restriction  operators  are  such  that 

pH  =  FHFh,  (4.24) 

then  PHU  =  PH(PhU)  and  (4.23)  reduces  to 

M{WH,FH(FhU)\H,c)-  M{WH,FhU]h,c)  =  B{WH ,PhU)  -  B(WH ,PH {PhU)), 

VWH  eVHyt  e]Tu  T2[.  (4.25) 

Equation  (4.25)  is  the  variational  counterpart  of  the  Germano  identity.  It  involves  only  the 
optimal  solution  FhU  in  the  space  Vh.  It  is  easy  to  verify  that  the  condition  (4.24)  holds  for 
several  restriction  operators  including  L2  or  Hl  projectors  and  interpolation  operators. 

At  this  stage  it  is  instructive  to  compare  this  equation  with  its  filtered  counterpart  (4.9). 

1.  In  the  filtered  case  only  the  non-linear  terms  contribute  to  the  right  hand  side  of  the 
Germano  identity.  In  the  variational  case  every  term  contributes. 

2.  In  the  filtered  case  it  is  assumed  that  the  Germano  identity  is  satisfied  pointwise,  thus 
when  the  residual  of  (4.9)  is  multiplied  by  any  weighting  function  and  integrated  over 
G,  the  result  is  zero.  In  contrast  to  this  (4.19)  holds  only  for  weighting  functions 
contained  in  VH ,  indicating  that  there  is  less  usable  information  in  the  variational 
Germano  identity.  This  is  direct  consequence  of  the  fact  that  fewer  constraints  are 
placed  on  the  model  term  in  the  variational  approach  than  in  the  filtered  approach. 
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3.  In  the  filtered  case  the  Germano  identity  is  a  tensor  relation.  In  the  variational  case  it 
is  a  weak  form  of  a  vector  equation.  It  this  aspect  the  variational  Germano  identity  is 
closer  to  the  vector  level  identity  studied  in  [45]. 

It  is  easy  to  see  that  the  differences  described  above  can  be  traced  back  to  differences  in 
the  definition  of  the  ideal  model  term  in  the  variational  and  filtered  cases. 

While  utilizing  this  identity  in  a  numerical  method  we  will  replace  P hU  with  the  numerical 
solution  Uh.  In  that  case  the  modeled  system  is  given  by  (4.17)  and 

M(Wh,  P ttUh\  H,  c )  -  M(WH,  Uh;  h,  c)  =  B{WH ,  Uh)  -  B(WH ,  PHUh), 

\/WH  €  V"  Vt  €]Ii,  T2[.  (4.26) 

The  equation  above  may  then  be  used  to  determine  the  parameters  c  in  the  model  term. 
This  modeled  system,  that  is  (4.17)  and  (4.26),  has  the  following  special  property:  If  (4-%4) 
holds,  then  it  permits  as  its  solution  the  optimal  solutions  P hU  and  PHU  on  the  spaces  Vh 
and  VH  respectively. 

Equation  (4.26)  represents  as  many  equations  as  the  dimension  of  VH .  Very  often  the 
number  of  parameters  to  be  determined  (that  is  the  number  of  components  of  c)  is  much 
smaller  than  this  number.  Motivated  by  what  is  done  in  the  traditional  dynamic  models,  we 
propose  the  following  methods  to  reduce  the  dimension  of  these  equations. 

Dissipation  method  This  method  is  motivated  by  the  approach  described  in  [15]  for  the 
filter-based  Germano  identity.  It  involves  setting  WH  —  ¥HUh  in  (4.26)  to  arrive  at 

M(PHUh,  P HUh\  H,  c)  -  M(PHUh,  Uh\  h,  c)  =  B(PHUh ,  Uh)  -  B(PH Uh,  PHUh), 

VtE}TuT2{.  (4.27) 

The  equation  above  represents  a  scalar  relation  that  may  be  used  to  determine  a  single 
parameter  in  the  model. 

Least  squares  method  This  method  is  based  on  the  approach  developed  in  [27]  for  the 
filter-based  Germano  identity.  In  this  case  we  set  WH  —  <pA  ,  A  =  1,  •  •  •  ,  dim(V/l),  to  arrive 
at 

M(4>A,¥HUh',  H,  c)  —  M((f>A,  Uh\  h,  c)  =  B(<f>A,Uh)  -  B(<t>A,PHUh), 

A  =  , dim(V/l) , Vt  e]Ti , T2 [(4.28) 

Thereafter  we  evaluate  the  c  that  minimizes  the  sum  of  the  square  of  the  residual  of  (4.28) 
for  each  A. 


4.5  Decay  of  Homogeneous  Isotropic  Turbulence 

In  this  section  we  apply  the  variational  Germano  identity  to  calculate  the  time-dependent 
viscosity  parameter  for  the  Smagorinsky  model  during  the  decay  of  homogeneous  isotropic 
turbulence. 
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The  strong  form  of  the  problem  is  given  by  (4.1)  &  (4.3),  where  Q  =]0,  L[3  and  its 
boundary  comprises  of  six  faces,  1^(0),  Tj(L),  j  =  1,2,3,  where 

Tj(c)  =  {x  G  dfl  |  Xj  =  c}  (4.29) 

It  is  assumed  the  boundary  conditions  are  periodic,  that  is 

u(x  +  Lej,t)  =  u(x,t ),  x  G  Tj(0),  t  g]Ti,T2[  (4.30) 

where  e,-  is  the  Cartesian  basis  vector  in  the  Xj  direction. 

The  equivalent  weak  form  is  given  by  (4.11).  We  solve  this  problem  using  a  Fourier- 
spectral  discretization.  The  space  of  weak  solutions  and  weighting  functions,  V,  is  given 
by 

v={v  I  V=  ■£  Vkc‘kx,  v_k  =  vk; 

(Lk/2ir)ez3 

J2  l*fc!2(1  +  ifcl2)  <°°;  lffcl2<°°}-  (4-31) 

(Lk/ 2tt)€Z3  ( Lk/2iT)€Z 3 

In  the  equation  above,  V ^  are  the  Fourier-coefficients  of  V  for  a  given  wave- vector  k,  the 
superscript  *  denotes  the  complex-conjugate  of  a  quantity  and  for  a  vector  s,  |s|  =  y/s*~s. 
The  terms  within  the  curly  brackets  imply  that  velocity  fields  contained  in  V  are  real- valued 
and  belong  to  FT1  (SI)  (the  energy  and  the  enstrophy  of  the  flow  are  finite  for  t  e]TuT2{), 
and  that  the  pressure  is  real- valued  and  square  integrable.  Note  that  the  periodic  boundary 
conditions  are  also  built  into  the  definition  of  V. 

The  modeled  system  of  equations  is  given  by  (4.17),  where  the  finite  dimensional  space 
is  given  by 

Vh  =  {Ve  V;  Vk  =  0,  |fc|  >  kh}  (4.32) 

and  the  model  term,  which  represents  the  Smagorinsky  eddy  viscosity  model  [30],  is  given 
by 

M(Wh ,  Uh;  h,  c)  =  2(Clh)2(Vswh,  [ Vsuh|' Vsuh).  (4.33) 

Note  that  in  (4.32),  \k\  —  Vk*  ■  k.  The  unknown  parameter  c\  in  the  Smagorinsky  model  is 
to  be  determined  using  the  variational  Germano  identity.  It  is  assumed  that  ci  is  independent 
of  the  spatial  coordinate  x ,  and  varies  only  in  time.  We  consider  the  use  of  the  dissipation 
method  (4.27)  and  the  least  squares  method  (4.28)  in  determining  this  parameter.  In  order  to 
utilize  these  relations  we  first  define  the  space  VH  and  the  restriction  operator  Following 
(4.32)  we  have 

VH  =  {V  eV]Vk  =  0,\k\>  kH},  (4.34) 

where  H  >  h  and  kH  <kh.  We  select  FH  to  be  the  L2  projector  of  V  into  VH .  Thus  FHV, 
where  V  €  V,  is  given  by 

( WH ,  V)  =  (WH,  V),  \/WH  G  VH.  (4.35) 
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Using  (4.34)  and  the  orthogonality  of  the  Fourier  modes  in  fl,  from  the  above  equation  we 
have 


y keikx  =  ¥Hy>  (4-36) 

\k\  <  kH 

where  ¥H  represents  a  sharp  cut-off  filter  in  the  wavenumber  space.  Note  that  while  the 
operation  of  FH  is  defined  only  on  elements  of  V,  ¥H  may  operate  on  any  scalar,  vector  or 
tensor  with  a  finite  L2  norm. 

We  are  also  interested  in  comparing  the  expressions  for  the  Smagorinsky  parameter 
obtained  using  the  variational  counterpart  of  the  Germano  identity  with  those  obtained 
from  the  filtered  version  of  this  identity.  In  order  to  be  consistent  in  this  comparison  we 
choose  the  filtering  operator  to  be  the  sharp  cut-off  filter  in  wavenumber  space. 


Dissipation  method  We  are  now  in  a  position  to  estimate  the  Smagorinsky  parameter. 
We  first  consider  the  dissipation  method  applied  to  the  variational  Germano  identity.  Using 
(4.33),  (4.12)  and  (4.36)  in  (4.27)  we  arrive  at 


1  {SH,uH  ®uH)  -  (SH,uh<g>uh) 
2H2(SH,\SH\SI1)-h?{SH,\Sh\Shy  J  11  2l’ 


where  u11  =  F Huh,  Sh  =  Vs uh  and  SH  =  S7suH.  In  deriving  the  above  relation  we  have 
made  use  of  the  L2  orthogonality  of  Fourier  modes  to  eliminate  all  the  bilinear  terms  in 
(4.27). 

The  dissipation  method  when  applied  to  filtered  form  of  the  Germano  identity  yields  (see 
for  example  [15]), 


1  (Sh,uH  ®uH)  -  {SH,uh®uh) 
2H*{Sh,\SH\SH)-h2{SH,\Sh\Shy  J  11  2 


(4.38) 


In  comparing  (4.37)  and  (4.38)  the  following  observations  may  be  made: 

1.  For  both  the  variational  and  the  filtered  approach  ((4.37)  and  (4.38),  respectively)  the 
contribution  from  the  linear  terms  on  the  right  hand  side  of  the  Germano  identity  to 
the  numerator  in  the  expression  for  ci,  is  zero.  In  the  filtered  case,  this  is  a  consequence 
of  the  Germano  identity  itself,  whereas  in  the  variational  case  this  is  due  to  the  special 
orthogonal  properties  of  the  basis  functions  and  the  restriction  operator. 

2.  The  first  term  in  the  numerator  and  the  denominator  of  the  variational  equation  (4.37) 
contains  SH  in  the  weighting  function  slot,  which  is  replaced  by  Sh  in  the  filtered 
equation  (4.38).  This  is  the  only  difference  between  these  two  expressions  and  it  is  a 
direct  consequence  of  the  fact  that  in  the  variational  case  constraints  are  imposed  on 
the  ideal  model  only  for  weighting  functions  contained  in  a  finite  dimensional  subspace, 
whereas  in  the  filtered  case  constraints  are  imposed  for  all  weighting  functions. 
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3.  The  fact  that  the  filtered  Germano  identity  represents  a  tensor  relation  and  the 
variational  Germano  identity  represents  a  vector  relation,  does  not  lead  to  any 
differences  in  the  final  expression  for  C\  when  using  the  dissipation  method.  This 
is  because  in  the  filtered  case,  when  using  the  dissipation  method,  the  equation  for 
ci  is  obtained  by  contracting  the  Germano  identity  with  the  rate  of  strain  tensor  and 
integrating  the  resulting  scalar  equation  over  Q.  Using  integration  by  parts  it  can 
be  shown  that  this  is  identical  to  contracting  the  divergence  of  the  filtered  Germano 
identity  with  the  velocity  field  and  integrating  the  resulting  scalar  equation  in  Q,  which 
in  turn  amounts  to  interpreting  the  Germano  identity  as  a  vector  relation. 


Least  squares  method  In  order  to  estimate  Ci  using  the  least  squares  method  applied  to 
the  variational  form  of  the  Germano  identity  we  set 

WH  =  eje-ik  x,  j  =  1,2,3, 4;  |fc|  <  kH ,  (4.39) 

in  (4.26)  and  then  use  (4.33),  (4.12),  (4.36)  and  the  L2  orthogonality  of  Fourier  modes  to 
arrive  at 

2c?FhV-A4  =  F  hV-AT,  V*€]Ti,T2[.  (4.40) 

where 

M  =  h2\Sh\Sh-H2\SH\SH,  (4.41) 

M  =  uh  <8>  uh  -  u11  <g>  uH ,  (4.42) 

and  as  before  uH  =  F Huh,  Sh  =  Vsuh  and  SH  =  Vsun .  Selecting  ci  such  that  sum  of  the 
squares  of  the  residual  of  (4.40)  is  minimized,  we  arrive  at 


Cl 


/l  (FHV  •A/' 

V  2(F*V  •  A4,F*V  •  M)’Vt  GJ  u  al‘ 


(4.43) 


The  least  squares  approach  when  applied  to  the  filtered  from  of  the  Germano  identity 
yields  (see  [27]) 


where 


Ci  = 


11  (A f,M) 
2  (M,M) 


,VtG]Ti,T2[, 


(4.44) 


M  =  ¥H(h2\Sh\Sh)-H2\SH\SH,  (4.45) 

M  =  FH(uh  <g>  uh)  -  uH  ®  uH.  (4.46) 

In  comparing  the  expressions  for  the  parameter  derived  using  the  variational  and  the 
filtered  versions  of  the  least  squares  approach,  that  is  (4.43)  and  (4.44),  the  following 
observations  may  be  made: 


52 


1.  In  this  case  also  for  both  the  variational  and  the  filtered  approaches  the  contribution 
from  the  linear  terms  to  the  Germano  identity  is  zero.  In  the  filtered  case,  this  is  a 
consequence  of  the  identity  itself,  whereas  in  the  variational  case  this  is  due  to  the 
properties  of  the  basis  functions  and  the  restriction  operator. 

2.  The  definitions  of  J\f  and  M.  for  the  variational  and  filtered  cases  are  similar.  However, 
in  the  variational  case  in  equation  (4.43),  these  tensors  are  operated  upon  by  the  sharp 
cut-off  filter  F^,  whereas  in  the  filtered  case  ((4.45)  &  (4.46)),  ¥H  acts  on  only  one 
term  in  these  tensors.  This  difference  may  be  traced  back  to  the  observation  that  in 
the  variational  approach  the  ideal  model  is  constrained  to  satisfy  a  given  relation  only 
for  a  finite  subset  of  weighting  functions,  whereas  in  the  filtered  case  it  is  required  to 
satisfy  the  relation  for  every  weighting  function. 

3.  The  inner  products  that  appear  in  the  definition  of  Ci  in  the  variational  case  (4.43) 
involve  the  divergence  of  the  tensors  H  and  M,  whereas  in  the  filtered  case  (4.44) 
they  involve  the  tensors  themselves.  This  is  because  the  expression  for  the  ideal  model 
and  the  Germano  identity  are  vector  relations  in  the  variational  case,  whereas  they  are 
tensor  relations  in  the  filtered  case. 

4.6  Numerical  Examples 

In  this  section  we  evaluate  the  performance  of  the  variational  Germano  identity  in  predicting 
the  time  dependent  Smagorinsky  constant  during  the  large  eddy  simulation  of  the  decay  of 
homogeneous  isotropic  turbulence.  We  compare  the  following  methods: 

1.  The  Dynamic  Smagorinsky  method,  where  the  Smagorinsky  parameter  ci  is  obtained 
using  the  dissipation  method  applied  to  the  variational  Germano  identity  (4.37). 

2.  The  Dynamic  Smagorinsky  method  with  the  least  squares  form  of  the  variational 
Germano  identity  (4.43). 

3.  The  Dynamic  Smagorinsky  method  with  the  dissipation  form  of  the  filtered  Germano 
identity  (4.38). 

4.  The  Dynamic  Smagorinsky  method  with  the  least  squares  form  of  the  filtered  Germano 
identity  (4.44). 

In  the  first  problem  we  consider  the  decay  of  turbulence  from  a  Taylor  micro-scale 
Reynolds  number  (Rex)  of  90  to  about  60  and  compare  the  LES  solutions  to  a  benchmark 
DNS  calculation.  In  the  second  problem  Rex  decays  from  716  to  626  and  we  compare  the 
LES  results  with  experimental  data.  In  both  cases  we  perform  the  simulations  using  a 
Fourier-spectral  code.  We  advance  in  time  using  a  third-order  Runge-Kutta  scheme.  The 
nonlinear  terms  (convective  acceleration  term  and  the  model  term)  are  evaluated  on  a  grid 
of  size  (31V/2)3.  The  effect  of  molecular  viscosity  is  accounted  for  exactly  by  utilizing  an 
appropriate  integrating  factor  for  each  ODE  [33]. 
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4.6.1  Rex  =  90 

We  compute  the  DNS  solution  with  2563  modes.  For  this  solution  we  use  a  random  initial 
condition  with  an  energy  spectrum  given  by  E(k)  =  (q2/2A)(kA/kp)  exp(— 2(k/kp)2),  where 

the  initial  turbulent  kinetic  energy,  2"  =  f ,  and  the  spectrum  is  peaked  at  kp  =  3.  The 
length  of  the  side  of  the  cube,  L  =  2ir  and  the  kinematic  viscosity  u  =  0.005.  The  same 
spectrum  was  used  in  [31]  to  study  the  decay  of  low  Reynolds  number  homogeneous  isotropic 
turbulence.  The  phase  of  the  modes  for  the  initial  velocity  field  is  chosen  randomly  and  each 
mode  satisfies  the  divergence-free  condition.  The  solution  is  allowed  to  evolve  according 
to  the  Navier-Stokes  equations  till  a  spectra  with  a  physical  A;-5/3  range  is  obtained  at 
t  =  T\  «  2.2.  The  evolution  of  the  kinetic  energy  on  a  cross-section  of  the  2563  DNS  is  shown 
in  the  video  clip.  The  animation  starts  from  the  random  initial  condition  and  ends  at  t  =  T\. 
At  this  instant  the  Taylor-microscale  Reynolds  number  is  «  90.  This  velocity  field  is  used 
as  an  initial  condition  for  all  LES  calculations,  which  are  performed  with  kh/(2n/L)  =  32 
and  kh/kH  =  2.  The  DNS  and  LES  solutions  are  evolved  till  t  =  T2  —  3  which  corresponds 
to  Re\  =  60  and  the  results  compared. 

In  Figure  4.1,  we  have  plotted  the  variation  of  the  Smagorinsky  parameter  c\  as  a  function 
of  time  for  the  LES  methods.  We  observe  that  coefficient  predicted  by  the  filtered  dissipation 
approach  is  much  higher  than  all  the  other  coefficients.  The  coefficient  predicted  by  the 
filtered  least  squares  approach  is  smaller  than  the  dissipation  approach.  We  also  observe  that 
the  coefficient  predicted  by  the  variational  dissipation  approach  is  very  close  to  the  filtered 
least  squares  approach  and  that  the  coefficients  predicted  by  the  two  variational  methods 
(dissipation  and  least  squares)  are  closer  to  each  other  when  compared  to  the  coefficients 
predicted  by  the  two  filtered  methods. 

It  is  worth  pointing  out  that  all  the  LES  methods  considered  here  differ  from  each  other 
only  in  the  way  the  Smagorinsky  parameter  is  evaluated.  Hence  the  differences  observed  in 
the  values  of  this  parameter  are  solely  responsible  for  any  difference  in  the  final  solutions. 

In  Figure  4.2,  we  have  plotted  the  variation  of  the  total  resolved  kinetic  energy  (|fc|  <  32) 
for  the  benchmark  DNS  solution  and  the  LES  solutions  as  a  function  of  time.  We  observe 
that  to  begin  with  all  the  models  are  overly  dissipative.  However  they  improve  with  time. 
We  also  observe  that  the  filtered  dissipation  approach,  which  attains  the  largest  viscosity 
(see  Figure  4.1),  is  the  most  dissipative,  whereas  the  performance  of  all  the  other  methods 
is  virtually  indistinguishable. 

In  Figure  4.3,  we  compare  the  energy  spectra  at  t  =  T2  =  3.  We  observe  that  all  the 
methods  under-predict  the  energy  in  the  fine  scale  modes.  This  in  turn,  lowers  the  transfer 
of  energy  from  the  coarse  scale  modes  to  the  fine  scale  modes  leading  to  an  over-prediction 
of  the  energy  in  the  coarse  scale  modes.  We  observe  that  the  method  with  the  largest 
value  of  the  Smagorinsky  parameter,  namely  the  filtered  dissipation  approach,  is  the  least 
accurate,  whereas  the  one  with  the  smallest  parameter,  the  filtered  least  squares  approach  is 
the  most  accurate.  The  variational  dissipation  approach  yields  results  that  are  very  similar 
to  the  filtered  least  squares  approach,  whereas  the  variational  least  squares  is  somewhat  more 
inaccurate  (though  not  as  inaccurate  as  the  filtered  dissipation  approach). 
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4.6.2  Rex  =  716 

In  this  section  we  compare  the  performance  of  the  LES  methods  at  a  higher  value  of 
Reynolds  number.  We  utilize  the  experiment  reported  in  [46]  as  the  benchmark  solution. 
In  our  simulations,  the  length  of  the  side  of  the  cube  L  =  5.12,  the  kinematic  viscosity 
u  =  1.5074  x  10-4  and  the  initial  spectra  is  given  by  equation  (6)  in  [46].  A  divergence- 
free  velocity  field  with  random  phases  is  used  to  initialize  all  LES  calculations  which  are 
performed  with  kh/(2Tr/L)  =  32  and  kh/kH  =  2.  The  LES  calculations  are  allowed  to  evolve 
till  t  —  0.1267,  when  correlations  between  the  velocity  coefficients  across  different  modes 
have  been  established  and  the  Smagorinsky  parameters  predicted  by  the  various  dynamic 
procedures  have  attained  a  steady  non-zero  value.  The  velocity  field  is  then  rescaled  so  that 
the  spectra  attains  its  original  value.  This  velocity  field  is  used  initial  condition  for  the  LES 
simulations.  Using  this  approach  the  effect  of  the  transient  associated  with  the  Smagorinsky 
parameter  starting  from  zero  (due  to  random  phases)  is  eliminated. 

In  Figure  4.4,  we  have  plotted  the  variation  of  the  Smagorinsky  coefficient  as  a  functidn 
of  time.  Note  that  in  this  figure  time  is  measured  in  the  x/M  units  used  in  [46].  We  note 
that  when  compared  with  Figure  4.1,  the  values  of  the  parameters  in  this  figure  are  generally 
higher.  This  may  be  attributed  to  a  more  extensive  inertial  range  for  the  simulation  reported 
in  Figure  4.4.  We  also  observe  that  as  before  the  filtered  dissipation  approach  predicts  the 
largest  value  of  the  Smagorinsky  parameter  and  the  filtered  least  squares  approach  predicts 
the  smallest  value.  The  variational  methods  predict  values  that  are  very  close  to  each  other. 
These  value  are  closer  to  the  value  for  filtered  least  squares  approach  than  to  the  filtered 
dissipation  approach. 

In  Figure  4.5,  we  have  plotted  the  Energy  spectra  at  t  —  T2  =  0.506  for  the  LES 
simulations  and  the  experimental  result.  This  time  corresponds  to  x/M  —  48  in  Kang 
et  al’s  experiment.  We  observe  that  the  filtered  dissipation  approach  has  significantly 
under-predicted  the  spectra  at  high  wavenumbers.  On  the  other  hand,  the  filtered  least 
squares  approach  is  somewhat  under-dissipative  leading  to  a  slight  pile-up  of  energy  at  high 
wavenumbers.  The  results  for  the  two  variational  methods  (dissipation  and  least  squares)  lie 
in  between  the  two  filtered  methods  and  are  about  as  accurate  as  the  filtered  least  squares 
method. 

4.6.3  Summary  and  explanation  of  numerical  results 

Based  on  the  results  presented  in  this  section  (at  Re \  —  90  and  716)  the  following  comments 
may  be  made  regarding  the  variational  and  the  filtered  forms  of  the  Germano  identity. 

1.  The  variational  form  yields  values  of  parameters  that  are  less  sensitive  to  the  method 
utilized  in  calculating  them  (dissipation  or  least  squares). 

2.  The  dissipation  approach  is  observed  to  yield  accurate  results  in  the  variational  case, 
whereas  in  the  filtered  case  it  is  overly  dissipative. 

3.  The  results  for  the  variational  dissipation  approach  are  about  as  accurate  as  that  of 
the  filtered  least  squares  approach.  However  the  costs  associated  with  implementing 
this  approach  are  lower.  Thus  it  may  represent  a  new,  easier  method  for  evaluating  a 
“more  accurate”  Smagorinsky  parameter. 
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From  the  numerical  results  presented  in  the  previous  section  we  observe  that  the  value  of 
the  Smagorinsky  parameter  calculated  using  the  filtered  dissipation  approach  is  the  highest. 
We  now  provide  an  explanation  for  this.  We  conduct  an  a-priori  analysis  using  the  DNS 
field  at  Re\  —  90.  We  compare  the  expression  for  the  filtered  dissipation  approach  with  its 
variational  counterpart.  The  difference  between  the  numerators  in  these  two  expressions  is 
given  by 

rkh 

-SH,uH  ®uH)=  /  Tn(k)dk  (4.47) 

Jk» 

where 

Tn{k)  =  -  [  u*k-  FT[V  •  ( uH  <g>  uH)]fcdrV  (4.48) 

J  r* 


(Sh 


and  k  =  |fc|.  In  the  integral  above,  T*  denotes  the  surface  of  a  sphere  of  radius  k  in  the 
wavenumber  space.  Similarly,  the  difference  between  the  denominator  in  (4.38)  and  (4.37) 
is  given  by 


where 


(Sh  -  SH,  |SH|SH)  =  I*  Td(k)dk 

JkH 

Uk)  =  -/**.-  ft[v  •  (|s''|s")]J.dri. 


(4.49) 
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In  Figure  4.6,  we  have  plotted  Tn  and  Td,  scaled  by  their  maximum  values,  as  a  function  of 
k.  We  observe  that  Td  rapidly  drops  to  zero  beyond  k11 .  As  a  result  the  denominators  in 
the  expressions  for  c\  for  the  filtered  and  the  variational  dissipation  approach  are  about  the 
same.  On  the  other  hand,  the  decay  in  T„  is  much  more  gradual  and  the  numerator  for  the 
filtered  approach  is  substantially  larger.  This  leads  to  a  larger  value  of  C\  for  the  filtered 
formulation. 

It  is  worth  noting  that  the  cause  of  this  discrepancy  between  the  filtered  and  the 
variational  formulations  is  the  difference  in  the  definition  of  an  ideal  model.  In  the  filtered 
formulation,  an  ideal  model  is  assumed  to  represent  the  subgrid  term  pointwise,  that  is  for 
all  weighting  functions.  This  leads  to  a  Germano  identity  wherein  contribution  from  a  model 
applied  at  a  cutoff  wavenumber  of  kH ,  to  modes  greater  than  kH ,  is  also  included.  This  is  not 
the  case  for  the  variational  formulation.  In  this  case  an  ideal  model  represents  the  effect  of 
the  subgrid  terms  only  on  a  select  number  of  weighting  functions  with  wavenumbers  smaller 
than  the  cutoff  wavenumber.  As  a  result  there  are  no  terms  in  the  Germano  identity  which 
contain  the  effect  of  the  model  beyond  the  cutoff  wavenumber.  We  note  that  this  difference 
between  these  two  approaches  is  present  for  all  flows  and  will  carry  over  to  more  complex, 
inhomogeneous  flows. 

For  the  least  squares  approach,  we  observe  that  the  difference  in  the  Smagorinsky 
parameter  evaluated  using  the  filtered  and  the  variational  formulations  is  smaller.  This 
is  remarkable  considering  that  the  terms  responsible  for  the  differences  observed  in  the 
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dissipation  approach,  are  also  present  in  the  least  squares  approach  ( uH  <g>  uH  in  the 
numerator  of  the  filtered  formulation  versus  F H(uH®uH)  in  the  numerator  of  the  variational 
formulation).  However  there  is  one  crucial  difference.  The  numerators  in  the  least  squares 
approach  contain  the  inner  product  of  these  terms  with  “other”  terms,  whose  Fourier 
coefficients  beyond  kH  have  been  verified  to  be  either  very  small  (for  the  filtered  formulation) 
or  identically  zero  (for  the  variational  formulation).  As  a  result,  in  this  case,  the  differences 
in  uH  ®  uH  and  F H{uH  <g>  uH)  beyond  kH  do  not  play  a  significant  role  determining  the 
Smagorinsky  parameter. 

It  is  important  to  bear  in  mind  that  in  the  least  squares  approach,  the  divergence  operator, 
which  is  present  in  the  variational  formulation  and  absent  from  the  dissipation  formulation 
is  another  source  of  differences.  The  precise  effect  of  this  operator  on  the  Smagorinsky 
parameter  will  be  investigated  in  the  future. 

4.7  Conclusions 

In  this  chapter  we  have  derived  the  variational  counterpart  of  the  Germano  identity  for 
the  incompressible  Navier-Stokes  equations.  This  identity  provides  a  means  of  determining 
model  parameters  within  the  variational  formulation  of  LES  and  is  analogous  to  the  Germano 
identity  derived  in  [15]  for  the  filter-based  LES  formulation.  We  have  identified  the  following 
differences  among  the  variational  and  the  filtered  versions  of  the  Germano  identity  and  traced 
their  origin  to  differences  in  the  definition  of  an  ideal  model  in  the  two  cases: 

1.  Linear  terms  in  the  Navier-Stokes  equations  contribute  to  the  variational  formulation 
and  not  to  the  filter-based  formulation. 

2.  The  variational  formulation  holds  only  for  a  finite  dimensional  space  of  weighting 
functions,  whereas  the  filter-based  formulation  holds  pointwise. 

3.  The  variational  formulation  is  essentially  a  vector  relation,  while  the  filter-based 
formulation  is  a  tensor  relation. 

In  light  of  the  third  point  above,  we  have  noted  that  the  variational  Germano  identity  is 
similar  to  (though  not  the  same  as)  the  vector- level  Germano  identity  discussed  in  [45]. 

We  have  applied  the  variational  Germano  identity  to  determine  the  Smagorinsky 
parameter  in  the  decay  of  homogeneous  isotropic  turbulence  and  compared  its  performance 
with  the  filter-based  formulation.  We  have  found  it  to  be  accurate  and  robust,  in  that  the 
value  of  the  parameter  is  less  sensitive  to  the  exact  approach  (dissipation  or  least-squares) 
used  in  evaluating  it.  The  example  described  in  this  chapter  represents  an  initial  result 
of  the  application  of  the  variational  Germano  identity.  Other  interesting  applications  to  be 
pursued  in  the  future  include  determining  multiscale  parameters  in  the  variational  multiscale 
formulation  and  determining  parameters  in  cases  (such  as  when  using  the  finite  element 
method)  where  contributions  from  linear  terms  are  non-zero. 
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Figure  4.1:  Variation  of  the  Smagorinsky  parameter  c\  with  time  for  the  low  Reynolds 
number  simulation. 
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Figure  4.2:  Variation  of  the  resolved  kinetic  energy  (\k\  <  32)  with  time  for  the  low  Reynolds 
number  simulation. 
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Figure  4.3:  Energy  spectra  at  t  =  T2  =  3  for  the  low  Reynolds  number  simulation. 
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Figure  4.4:  Variation  of  the  Smagorinsky  parameter  ci  with  time  for  the  high  Reynolds 
number  simulation. 


Chapter  5 


Variational  Germano  Identity  Applied 
to  Computing  Hydrodynamic  Noise 

5.1  Introduction 

In  this  chapter  we  assess  the  performance  of  the  LES  models  developed  as  part  of  this  research 
program  in  computing  far-field  noise.  We  consider  the  simple  case  of  homogeneous  isotropic 
turbulence  and  apply  Lighthill’s  acoustic  analogy  to  compute  noise.  This  decouples  the 
problem  into  two  parts.  In  the  first  part  we  solve  the  incompressible  Navier  Stokes  equations 
to  determine  the  fluctuating  (turbulent)  velocity  field.  For  this  purpose  we  consider  sustained 
turbulence  which  is  driven  by  a  forcing  at  low  wavenumbers.  We  solve  this  problem  on  a 
fine  mesh,  where  the  dissipation  length  scale  is  resolved,  and  then  on  a  much  coarser  mesh 
where  various  LES  models  are  employed  to  represent  the  effect  of  the  missing  scales. 

The  velocity  field  from  the  turbulent  problem  is  used  to  construct  quadrupole  sources 
which  drive  the  acoustic  problem.  The  acoustic  problem  is  transformed  into  the  frequency 
domain  and  then  solved.  Since  we  are  dealing  with  an  unbounded  domain  with  simple 
sources,  an  analytical  solution  to  this  problem  is  feasible.  This  solution  is  constructed 
and  used  to  determine  the  far-field  acoustic  pressure  intensity  as  a  function  frequency  and 
wavenumber.  The  acoustic  intensity  is  computed  for  sources  obtained  from  the  highly 
resolved  turbulent  simulation  (the  Direct  Numerical  Solution,  DNS)  and  those  obtained  from 
LES  models.  The  acoustic  intensity  for  the  LES  models  is  compared  to  that  of  the  DNS  and 
conclusions  about  their  efficacy  are  drawn.  This  rather  simple  framework  provides  a  useful 
means  to  categorize  the  performance  of  LES  models  in  predicting  turbulence  generated  noise. 

In  the  following  section  we  derive  an  expression  for  the  far-field  acoustic  intensity 
corresponding  to  a  homogeneous  isotropic  turbulence  field.  Thereafter  we  present  results 
comparing  the  performance  of  various  LES  models. 

5.2  Expression  for  the  Far-Field  Acoustic  Intensity 

The  use  of  LES  velocity  data  in  Light  hill’s  acoustic  analogy  allows  for  the  computation  of 
the  far-field  acoustic  pressure  for  low  Mach  number  flows.  The  effect  of  the  hydrodynamic 


fluid  motion  appears  as  a  quadrapole  source  term  in  the  wave  equation 

(^-y2)p  -  V'V'T 

t  =  u<S>u  (5.2) 

where  c  is  the  speed  of  sound  in  the  ambient  medium.  The  Helmholtz  equation  is 
obtained  by  taking  the  Fourier  transform  in  time  of  the  wave  equation,  such  that  p(x,  t )  = 
fZoPMe-^dw. 

-(fc2  +  V2)p  =  V-V-f  (5.3) 

where  f  is  the  Fourier  transform  in  time  of  r  and  k  =  -. 

C 

Writing  the  solution  in  terms  of  a  Green’s  function,  the  pressure  becomes: 


p(x,  u>)=  G(x;  y)Vy  ■  Vy  ■  f(y,  u;)dQv 

J 


(5.4) 


For  an  unbounded  domain,  the  Green’s  function  is 

1  e-*fc|x-i/| 


G(x\  y)  = 


An r  \x  —  y\ 


(5.5) 


The  integral  in  (5.4)  can  be  approximated  assuming  that  a:  is  in  the  far-field.  The  acoustic 
pressure  as  a  function  of  A;,  a;  and  x  is  then: 


p  «fc|x| 

p(x,  k,un)  =  -k2-r~YT^  ®  £  '■  T(k£,  u>)  (5.6) 

47T|a:| 

where  £  =  ^  is  the  unit  vector  parallel  to  x  and  T(/c£,  to)  is  the  space-time  Fourier  transform, 
with  wavenumber  /c£  and  frequency  u,  of  r  =  u®u  calculated  from  the  incompressible  velocity 
field.  This  quantity  is  readily  available  from  the  incompressible  Navier-Stokes  calculation 
provided  a  long  enough  time  series  of  velocity  data  can  be  stored  for  calculating  the  Fourier 
transform  in  time  of  the  source  term. 

Since  the  turbulence  is  isotropic,  an  ensemble  average  of  the  acoustic  pressure  is  equivalent 
to  an  average  over  all  directions,  £,  of  one  realization  of  the  acoustic  pressure.  By  averaging 
the  acoustic  pressure  over  all  directions,  a  two  variable  function  of  the  scaled  noise  intensity 
as  a  function  of  k  and  u>  is  obtained. 


(4n\x\)2  (p2(x,k,u))  =  (((A£)®(fc£)  :  T(k^UJ)f)  (5-7) 

That  is  the  the  scaled,  far-field  acoustic  intensity  is  determined  by  the  quantity 

(m)®m--T(kz,u;))2). 
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loq<P{k,o))2 


5.3  Evaluation  of  LES  Models 

In  Figure  5.1,  a  comparison  is  made  between  DNS  and  several  LES  schemes  with  regards 
to  the  scaled  acoustic  intensity.  Figure  5.1(a)  shows  the  results  from  the  DNS  velocity 
field  while  5.1(d)  shows  the  result  from  a  LES  with  a  constant  Smagorinsky  parameter.  A 
large  improvement  is  observed  when  the  variational  Germano  identity  is  used  to  determine 
the  Smagorinsky  parameter  (see  Figure  5.1(b)).  The  best  results  are  obtained  when 
the  variational  Germano  identity  is  used  to  determine  the  parameter  of  a  one  parameter 
multiscale  model  where  the  model  is  only  applied  to  fine  scale  Fourier  modes  (Figure  5.1(c)). 


log<P{k,oo)2>  -  DNS 


(a)  DNS 


■>  -  variational  all-small 


(b)  variational  dynamic  Smagorinsky 


k>g<P(k,!»)2>  -  variational  all-all 


(c)  variational  multiscale 


(d)  static  Smagorinsky 


Figure  5.1:  Comparison  between  DNS  and  LES  of  the  normalized  acoustic  pressure  intensity 
as  a  function  of  \k\  and  u. 


From  the  plots  in  Figure  5.1,  the  acoustic  pressure  frequency  spectrum  for  a  particular 
medium  is  obtained  by  choosing  a  sound  speed  and  taking  a  section  of  the  plot  along  k  =  u/c. 
In  each  plot,  lines  are  added  corresponding  to  values  of  c  =  10  and  50.  These  values  of  c 
correspond  to  a  Mach  number  of  about  M  —  0.1  and  M  =  0.02  which  is  well  within  the 


low-Mach  number  limit. 

When  an  LES  is  performed,  the  Lighthill  source  term  is  only  partially  known.  The 
contribution  from  the  fine  scales  is  not  available.  To  correct  for  this,  the  modeled  subgrid 
stress  has  been  included  in  the  source  term  such  that,  for  the  Smagorinsky  model,  r  = 
(uh  ®uh  +  2V2(C3h)2\Vsuh\Vsuh).  In  the  results  described  above  we  have  included  this 
term. 
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